plus sc10x
load dependancies
load 10x data
GEX.seur <- readRDS("./GEX0811.seur.preAnno.rds")
GEX.seur## An object of class Seurat
## 18371 features across 25338 samples within 1 assay
## Active assay: RNA (18371 features, 1229 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapcolor.FB <- ggsci::pal_igv("default")(49)[c(2,13,33,1,15,
34,26,28,40,18)]
color.cnt <- color.FB[c(10,9,7,
5,4,1)]DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno2")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "preAnno1")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "preAnno1", split.by = "DoubletFinder0.05")## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "FB.info", split.by = "DoubletFinder0.05")
advanced filtering
only keep DF0.05 Singlets
remove C1
GEX.seur <- subset(GEX.seur, subset = DoubletFinder0.05=="Singlet" & seurat_clusters %in% setdiff(c(0:19),c(15,18,19)))
GEX.seur <- subset(GEX.seur, subset = nFeature_RNA < 3000 & nCount_RNA < 8000)
GEX.seur## An object of class Seurat
## 18371 features across 23692 samples within 1 assay
## Active assay: RNA (18371 features, 1229 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapVlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "preAnno1")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "preAnno1")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.1, group.by = "FB.info")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "FB.info")
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "cnt", cols = color.cnt)
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "cnt", cols = color.cnt)
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "cnt", cols = color.cnt)
plota + plotb + plotc
GEX.seur$FB.new <- factor(as.character(GEX.seur$FB.info),
levels = setdiff(levels(GEX.seur$FB.info),
c("Doublet","Negative"))[c(10,8:9,6:7,
5,3:4,1:2)])
color.FBnew <- color.FB[c(10,8:9,6:7,5,3:4,1:2)]par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.new)
barplot(sl_stat,ylim = c(0,4200),
#col = c("#FF6C91","lightgrey",color.FB),
col = c(color.FBnew),
main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:10*1.2-0.48 ,levels(GEX.seur@meta.data$FB.new), las=3, cex.axis=0.85)
text(x=1:10*1.2-0.45,y=sl_stat+285,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
new
re-clustering
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1800)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 300)## [1] "Spp1" "Ccl4" "Cxcl10" "Ccl5"
## [5] "Ccl3" "Spp1-EGFP" "Ccl12" "Hist1h1b"
## [9] "Mki67" "Cxcl9" "Ifi27l2a" "Ccl2"
## [13] "Cxcl13" "Apoe" "Rsad2" "Top2a"
## [17] "Xist" "Fn1" "Cst7" "Hist1h2ae"
## [21] "Cenpf" "Lpl" "Ifitm3" "Hmgb2"
## [25] "Fabp5" "Hist1h2ap" "Lgals3" "Ifit3"
## [29] "Ch25h" "Ifit2" "Stmn1" "Prc1"
## [33] "Iigp1" "Isg15" "Ifit1" "Pclaf"
## [37] "Ccl7" "Cd52" "Hmmr" "Ube2c"
## [41] "Egr1" "Gpnmb" "Nusap1" "Mmp12"
## [45] "Cenpa" "Cd72" "Cd69" "Il1b"
## [49] "Apoc1" "Csf1" "Gm26737" "Aspm"
## [53] "Cd9" "Ctsd" "Cenpe" "Birc5"
## [57] "Cd74" "H2-K1" "Cdca8" "Pbld1"
## [61] "Ifi204" "Gdf15" "Lgals3bp" "H2-D1"
## [65] "H2afz" "Igf1" "Lyz2" "Hist1h3c"
## [69] "Tpx2" "Ccl6" "Slpi" "Ctsb"
## [73] "Saa3" "B230312C02Rik" "Arg1" "Il12b"
## [77] "Fth1" "Bst2" "Serpina3g" "Smc2"
## [81] "Mif" "Kif11" "Lgals1" "Acod1"
## [85] "Ccl8" "Esco2" "Fxyd5" "Ttr"
## [89] "B2m" "Smc4" "H2afx" "Knl1"
## [93] "Hist1h1e" "Ftl1" "Ccnb1" "Cd5l"
## [97] "Cdkn1a" "Ms4a7" "Cd63" "Pbk"
## [101] "Rrm2" "Mrc1" "Hmox1" "Igfbp5"
## [105] "Ifi211" "Slfn5" "Cybb" "Hba-a2"
## [109] "Hba-a1" "Kif15" "Gm26885" "Dqx1"
## [113] "Rcan1" "Hist2h2ac" "Wfdc17" "Kif23"
## [117] "Ctsz" "Mis18bp1" "Lilrb4a" "Ifi207"
## [121] "Cdca3" "Kif20b" "Ccnb2" "Apoc2"
## [125] "Mt1" "4933407L21Rik" "3110039M20Rik" "Gbp2"
## [129] "Syngr1" "H2-Q7" "Atad2" "Ifi203"
## [133] "Ctsl" "Ccr1" "Enpp2" "Marco"
## [137] "Gm15987" "Pf4" "Gm42047" "Cmpk2"
## [141] "Cdc20" "Ifi44" "Ifit3b" "Bub1b"
## [145] "Ccl9" "Cxcl2" "Tnfsf18" "Anln"
## [149] "Racgap1" "Dennd5b" "Aldh1a3" "Gfod2"
## [153] "Wdcp" "Atf3" "Sgo2a" "Tspo"
## [157] "Ifitm1" "Ckap2l" "Tc2n" "Cspg4"
## [161] "Spc24" "Ndc80" "Gm43409" "Plac8"
## [165] "Irf7" "Tubb4b" "Oasl2" "Ctla2a"
## [169] "Rbms3" "Myh4" "Lmnb1" "C4b"
## [173] "Gnas" "H2-Ab1" "Hist1h4d" "2210409D07Rik"
## [177] "Sh3tc2" "Cdk1" "Cox6a2" "Gapdh"
## [181] "Gm12863" "Dcx" "Timp2" "Kcnq1ot1"
## [185] "Phf11b" "Spc25" "Ly6a" "Ccl24"
## [189] "Ank" "Hist1h2ab" "Cdca2" "Ly86"
## [193] "Baiap2l2" "Fos" "Tk1" "Tsix"
## [197] "Lig1" "Cks2" "Cstb" "Ifi213"
## [201] "Clspn" "Ccnf" "Nfasc" "Gfap"
## [205] "Lockd" "Cks1b" "Cp" "Snca"
## [209] "Gm4951" "Cxcl14" "Tex11" "Fam20c"
## [213] "Ckap2" "Parp14" "Cxcl16" "Lrrc31"
## [217] "Gm10457" "Tmpo" "Ccna2" "Ifi209"
## [221] "Plk1" "Rbpms2" "Lsamp" "Ifi205"
## [225] "Loxl2" "Vim" "Tnr" "Prokr1"
## [229] "Ifi208" "Apob" "Ifnb1" "Fbxo5"
## [233] "Ncapg2" "Tlr2" "Incenp" "Pcna"
## [237] "Gm8251" "mt-Nd1" "Diaph3" "Ncapg"
## [241] "Zbp1" "Capg" "Emp3" "Ndnf"
## [245] "Gm35558" "Rnd3" "Rad51ap1" "F830016B08Rik"
## [249] "Apoc4" "Mamdc2" "Tubb5" "Aurkb"
## [253] "Trim59" "Crybb1" "Gadd45b" "Srgn"
## [257] "Clec4d" "Olfr433" "Fgl2" "Egr3"
## [261] "Gas2l3" "Kif4" "Nr4a3" "Il1rn"
## [265] "Slfn1" "Sgo1" "Stat1" "Gm43568"
## [269] "Cd83" "Oasl1" "Ifi206" "Gpr84"
## [273] "Tuba1b" "Usp18" "Id2" "Dtl"
## [277] "Dlgap5" "Mt2" "Arl6ip1" "Anxa5"
## [281] "H2-Aa" "mt-Co1" "Hist1h3e" "Ccdc34"
## [285] "Rgs17" "Olfr1369-ps1" "Rnase2a" "Cpd"
## [289] "E2f8" "Syt1" "Chp2" "Fam57b"
## [293] "Ctnna3" "Cbx5" "Rrm1" "Gm29773"
## [297] "Gdf3" "Nuf2" "Herc6" "Ier3"# exclude MT genes and more
# add sex-related Xist/Tsix
DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|Egr1|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AC|^AI|^AA|^AW|^AY|^BC|^Gm|^Hist|Rik$|-ps|Xist|Tsix|^Ifi|^Isg|^Mcm",
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)
VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
c(MT_gene,
DIG,
CC_gene) )
GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur), seed.use = 868)## PC_ 1
## Positive: Cst7, Cd72, B2m, Ch25h, Lgals3bp, Cd52, Fth1, Ctsb, Ctsz, Ccl3
## Cd9, Bst2, Fxyd5, Ctsd, H2-K1, Mif, Lpl, H2-D1, Fabp5, Spp1-EGFP
## Lilrb4a, Csf1, Tspo, Ccl2, Ccl4, Cd63, Lyz2, Zbp1, Oasl2, Fn1
## Negative: Crybb1, Ccr5, Sox4, Fchsd2, Hpgd, Clec4a3, Bank1, Fcgrt, Lst1, Klf7
## Mrc1, Igfbp4, Serpinf1, Camk2n1, Rgs7bp, Fry, Gpr165, Cbfa2t3, Tnfrsf17, Il7r
## Cryba4, Gbp7, Filip1l, Lrba, Ppm1l, Rdx, Zbtb20, Ccr1, Ddit4, Nav2
## PC_ 2
## Positive: Prc1, Pclaf, Pbk, Knl1, Esco2, Kif15, Aspm, Spc24, Stmn1, Neil3
## Smc2, Racgap1, Mis18bp1, Lockd, Bub1b, Ccna2, Sgo2a, Ccnb1, Ccnf, H2afx
## Plk1, Tk1, Cit, Spc25, Ncapg, Diaph3, Shcbp1, Trim59, Kif4, Fbxo5
## Negative: Ctsd, Ctsb, Tyrobp, B2m, H2-K1, Cst7, Cd9, Timp2, Fth1, H2-D1
## Cd63, Ctsz, Apoe, Lgals3bp, Lyz2, Ctsl, Ftl1, Grn, Npc2, Mpeg1
## Slfn5, Ly86, Ccl3, H2-Q7, H2-T23, Bst2, Oasl2, Syngr1, Cd72, Ccl4
## PC_ 3
## Positive: Rsad2, Iigp1, Slfn5, Oasl2, Parp14, Stat1, Usp18, Fgl2, Phf11d, Cxcl10
## Rtp4, Herc6, Slfn8, Trim30a, Ly6e, Ly6a, Slfn2, 9930111J21Rik2, Sp100, Oasl1
## Irf7, Rnf213, Zup1, Gbp2, Stat2, Tor3a, Xaf1, Ddx58, Ccl12, Etnk1
## Negative: Tmsb4x, Ftl1, Lpl, Tyrobp, Igf1, Mt1, Fabp5, Nme2, Serf2, Mif
## Rbm3, Myl6, Hint1, Cstb, Nme1, Spp1-EGFP, Gapdh, Lgals3, Cox6c, Atox1
## Emp3, Gpx1, Atp5g1, Gng5, Pkm, Gas6, Scd2, Calr, Fxyd5, Aldoa
## PC_ 4
## Positive: Ctsd, Cd9, Pmp22, Nav2, Samsn1, Fam102b, Hpgd, Csmd3, Lrba, Sgk1
## Ctsz, Cadm1, Tent5c, Syngr1, Thrsp, Pdgfb, Srgn, Lag3, Ccl6, Gpr183
## Cd34, Il7r, Ccl9, Prc1, Cd14, Ccnb1, Aspm, Timp2, Fth1, Nav3
## Negative: Tmsb4x, Mrc1, Rbm3, Rsad2, Iigp1, Ccl12, Igfbp4, Slfn5, Crybb1, Cxcl10
## Fcgrt, Gas6, Gbp7, Nme2, Cryba4, Oasl1, Oasl2, Usp18, Chd4, Gbp2
## Eif3a, Parp14, Ptma, Ccnd1, Lst1, Eif2ak2, Tpr, Hint1, Zbp1, Zfas1
## PC_ 5
## Positive: Tmsb4x, Apoe, Ccnb1, Aspm, Lyz2, Knstrn, Tyrobp, Cdkn3, Ftl1, Plk1
## Cep55, Ctsl, Zfas1, Cd52, Racgap1, Cox8a, Fth1, Apoc1, Lockd, Gpx1
## Prc1, H2afx, Bub1b, Tuba1c, B2m, Selenoh, Ramp1, Kif20a, Cox6c, Hacd4
## Negative: Lig1, Dhfr, Dut, Macf1, Atad5, Fignl1, Ccl4, Cdt1, Fam102b, Nsd2
## Dnmt1, Nav3, Wdhd1, Rif1, Kcnq1ot1, Rad54l, Rmi2, Topbp1, Timeless, Chaf1a
## Brca2, Ncapg2, Pmp22, Ccnd1, Tcf19, Pdgfb, Brca1, Nop56, Hnrnpab, Lilrb4alength(VariableFeatures(GEX.seur))## [1] 1389head(VariableFeatures(GEX.seur),300)## [1] "Spp1" "Ccl4" "Cxcl10" "Ccl5" "Ccl3" "Spp1-EGFP"
## [7] "Ccl12" "Cxcl9" "Ccl2" "Cxcl13" "Apoe" "Rsad2"
## [13] "Fn1" "Cst7" "Lpl" "Fabp5" "Lgals3" "Ch25h"
## [19] "Stmn1" "Prc1" "Iigp1" "Pclaf" "Ccl7" "Cd52"
## [25] "Gpnmb" "Mmp12" "Cd72" "Cd69" "Il1b" "Apoc1"
## [31] "Csf1" "Aspm" "Cd9" "Ctsd" "Cd74" "H2-K1"
## [37] "Pbld1" "Gdf15" "Lgals3bp" "H2-D1" "H2afz" "Igf1"
## [43] "Lyz2" "Ccl6" "Slpi" "Ctsb" "Saa3" "Arg1"
## [49] "Il12b" "Fth1" "Bst2" "Serpina3g" "Smc2" "Mif"
## [55] "Lgals1" "Acod1" "Ccl8" "Esco2" "Fxyd5" "Ttr"
## [61] "B2m" "H2afx" "Knl1" "Ftl1" "Ccnb1" "Cd5l"
## [67] "Cdkn1a" "Ms4a7" "Cd63" "Pbk" "Mrc1" "Hmox1"
## [73] "Igfbp5" "Slfn5" "Cybb" "Kif15" "Dqx1" "Rcan1"
## [79] "Wfdc17" "Ctsz" "Mis18bp1" "Lilrb4a" "Apoc2" "Mt1"
## [85] "Gbp2" "Syngr1" "H2-Q7" "Ctsl" "Ccr1" "Enpp2"
## [91] "Marco" "Pf4" "Cmpk2" "Bub1b" "Ccl9" "Cxcl2"
## [97] "Tnfsf18" "Racgap1" "Dennd5b" "Aldh1a3" "Gfod2" "Wdcp"
## [103] "Atf3" "Sgo2a" "Tspo" "Tc2n" "Cspg4" "Spc24"
## [109] "Plac8" "Irf7" "Oasl2" "Ctla2a" "Rbms3" "Myh4"
## [115] "Lmnb1" "C4b" "Gnas" "H2-Ab1" "Sh3tc2" "Cox6a2"
## [121] "Gapdh" "Dcx" "Timp2" "Kcnq1ot1" "Phf11b" "Spc25"
## [127] "Ly6a" "Ccl24" "Ank" "Ly86" "Baiap2l2" "Tk1"
## [133] "Lig1" "Cstb" "Ccnf" "Nfasc" "Gfap" "Lockd"
## [139] "Cp" "Snca" "Cxcl14" "Tex11" "Fam20c" "Parp14"
## [145] "Cxcl16" "Lrrc31" "Ccna2" "Plk1" "Rbpms2" "Lsamp"
## [151] "Loxl2" "Vim" "Tnr" "Prokr1" "Apob" "Ifnb1"
## [157] "Fbxo5" "Ncapg2" "Tlr2" "Incenp" "Diaph3" "Ncapg"
## [163] "Zbp1" "Capg" "Emp3" "Ndnf" "Rnd3" "Apoc4"
## [169] "Mamdc2" "Tubb5" "Trim59" "Crybb1" "Gadd45b" "Srgn"
## [175] "Clec4d" "Olfr433" "Fgl2" "Egr3" "Kif4" "Nr4a3"
## [181] "Il1rn" "Slfn1" "Sgo1" "Stat1" "Cd83" "Oasl1"
## [187] "Gpr84" "Tuba1b" "Usp18" "Id2" "Mt2" "Arl6ip1"
## [193] "Anxa5" "H2-Aa" "Ccdc34" "Rgs17" "Rnase2a" "Cpd"
## [199] "Syt1" "Chp2" "Fam57b" "Ctnna3" "Gdf3" "Herc6"
## [205] "Ier3" "Phf11d" "Cdk15" "Zdbf2" "Otor" "Cldn11"
## [211] "Sox2ot" "Lpar1" "Akr1b7" "Tril" "Fancd2os" "Fgf6"
## [217] "Vmn2r30" "Osbpl5" "Fat2" "Asic2" "Cyb561" "Sox9"
## [223] "F2rl1" "Blk" "Igsf11" "Ripply3" "Syt4" "Crtac1"
## [229] "Col4a5" "Cdkn3" "Lhfpl4" "Hif1a" "Knstrn" "Fabp3"
## [235] "H2afv" "Clmn" "Cst6" "Cck" "Ldha" "Neil3"
## [241] "Rnf213" "Slfn2" "Selenow" "Itgax" "Tpi1" "Ptger4"
## [247] "Alas2" "Il4i1" "Fcgr4" "Stil" "Ezh2" "Ms4a4c"
## [253] "Klhl33" "Tuba1c" "Pkm" "Hirip3" "Ecrg4" "Igfbp2"
## [259] "Il12rb1" "Cd24a" "Cd40" "Rgs1" "Olfr889" "Dut"
## [265] "Rab7b" "Rtp4" "Kmo" "Mcf2l" "Olr1" "Meg3"
## [271] "E2f7" "Anp32b" "Folr2" "Troap" "Adam33" "Cit"
## [277] "Gbp5" "Ly6e" "Shcbp1" "Melk" "Plau" "Zfp9"
## [283] "Oas3" "Colec12" "Cdkn2c" "Sag" "Otos" "Adcy10"
## [289] "Ermn" "Pde1a" "Grin2b" "Hs3st4" "Tex15" "Rasd2"
## [295] "Tox3" "Ccdc151" "Ptprk" "Ppp1r27" "Dpf3" "Sycp2l"DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.new", cols = color.FBnew, ) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.new", cols = color.FBnew)
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "cnt", cols = color.cnt)
cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2),
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)
cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1-EGFP", reduction = "pca",dims = 1:2),
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.new", split.by = "FB.new", cols = color.FBnew, ncol = 5)
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2, split.by = "cnt")
FeaturePlot(GEX.seur, features = "Spp1-EGFP", reduction = "pca",dims = 1:2, split.by = "cnt")
FeaturePlot(GEX.seur, features = "P2ry12", reduction = "pca",dims = 1:2, split.by = "cnt")
FeaturePlot(GEX.seur, features = "Mki67", reduction = "pca",dims = 1:2, split.by = "cnt")
cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "Spp1", reduction = "pca",dims = 1:2, split.by = "cnt", ) ,#& xlim(c(-7.5,33)) & ylim(c(-11,49)) ,
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "Spp1", reduction = "pca",dims = 1:2, split.by = "cnt") ,#& xlim(c(-7.5,33)) & ylim(c(-11,49)),
ncol = 1)
cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "Spp1-EGFP", reduction = "pca",dims = 1:2, split.by = "cnt") ,
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "Spp1-EGFP", reduction = "pca",dims = 1:2, split.by = "cnt") ,
ncol = 1)
cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "P2ry12", reduction = "pca",dims = 1:2),
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)
cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "P2ry12", reduction = "pca",dims = 1:2, split.by = "cnt") ,
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "P2ry12", reduction = "pca",dims = 1:2, split.by = "cnt") ,
ncol = 1)
cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Mki67", reduction = "pca",dims = 1:2),
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "cnt", cols = color.cnt),
rel_widths = c(5,5.6),
ncol = 2)
cowplot::plot_grid(
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P05.CTR","P05.MIG","P05.GFP")), features = "Mki67", reduction = "pca",dims = 1:2, split.by = "cnt") ,
FeaturePlot(subset(GEX.seur,subset=cnt %in% c("P28.CTR","P28.MIG","P28.GFP")), features = "Mki67", reduction = "pca",dims = 1:2, split.by = "cnt") ,
ncol = 1)
DimHeatmap(GEX.seur, dims = 1:12, cells = 6000, balanced = TRUE,ncol = 4,nfeatures = 80)
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(desc(PC_1) ) %>% head(40)## PC_1 PC_2 PC_3 PC_4 PC_5
## Cst7 0.13033723 -0.04914273 -0.024344645 0.0391614600 0.022557581
## Cd72 0.10462524 -0.03082428 -0.025466799 -0.0060531634 -0.020052539
## B2m 0.10373478 -0.05260759 0.042947295 0.0302277177 0.060053627
## Ch25h 0.10215044 -0.02902052 -0.016884600 -0.0136793196 -0.058608146
## Lgals3bp 0.10039744 -0.04596373 0.067935444 0.0083731503 0.024409626
## Cd52 0.10030771 -0.03034810 -0.012582220 -0.0374747502 0.072467231
## Fth1 0.10017184 -0.04848443 -0.014527262 0.0583658027 0.070162498
## Ctsb 0.09965710 -0.05693450 -0.041613614 0.0405126748 0.037755196
## Ctsz 0.09570110 -0.04719809 -0.016882643 0.0814963494 0.019413379
## Ccl3 0.09494543 -0.03390996 -0.033316477 0.0319485657 -0.066979156
## Cd9 0.09446644 -0.04879779 -0.007663390 0.1393953163 -0.020040494
## Bst2 0.09385601 -0.03219153 0.043282808 -0.0558837565 0.045345922
## Fxyd5 0.09371542 -0.01947810 -0.057627566 0.0095793204 0.004101436
## Ctsd 0.09312216 -0.07980194 0.010063702 0.1698541195 -0.012359922
## H2-K1 0.09307069 -0.04946374 0.061854974 0.0557680327 0.040462048
## Mif 0.09164890 -0.01956187 -0.072324527 0.0098647201 -0.003255225
## Lpl 0.08980026 -0.01488866 -0.090240695 -0.0003338312 -0.040050573
## H2-D1 0.08873270 -0.04841797 0.060894356 0.0517284292 0.055801786
## Fabp5 0.08762805 -0.02600469 -0.079178279 0.0150022563 -0.010474565
## Spp1-EGFP 0.08739725 -0.02331045 -0.066788887 0.0056088220 -0.017086488
## Lilrb4a 0.08727945 -0.02537670 -0.048874340 0.0171707711 -0.067661249
## Csf1 0.08626413 -0.02907911 -0.024130726 0.0133250495 -0.058306480
## Tspo 0.08615673 -0.02259411 0.020394771 -0.0420595092 0.055264464
## Ccl2 0.08516958 -0.02556872 0.044195824 -0.0586610844 -0.033474700
## Ccl4 0.08508932 -0.03059175 -0.024160255 0.0181325362 -0.085269328
## Cd63 0.08310536 -0.04724694 -0.054186027 0.0300100818 0.030584581
## Lyz2 0.08260736 -0.03988305 -0.022641236 0.0189047277 0.089614900
## Zbp1 0.08228250 -0.02797594 0.070085144 -0.0694688088 0.013992002
## Oasl2 0.08192573 -0.03214081 0.116312337 -0.0799565678 0.009603802
## Fn1 0.08104928 -0.02006132 -0.040377111 -0.0127196754 -0.013037203
## Irf7 0.08059542 -0.02868662 0.086758150 -0.0685628905 0.015186291
## Apoe 0.07978716 -0.04636109 -0.025883130 -0.0538887914 0.137889836
## H2-Q7 0.07970411 -0.03245787 0.052403842 0.0215995009 0.030129185
## Rsad2 0.07945714 -0.02853418 0.127246449 -0.1123588025 -0.011752638
## Spp1 0.07916311 -0.02501669 -0.056085597 0.0140606130 -0.004654707
## Cxcl16 0.07808397 -0.02833363 -0.001287169 0.0348527122 -0.019926816
## Phf11b 0.07709413 -0.02253058 0.070077143 -0.0586345152 0.009170246
## Ftl1 0.07580178 -0.03755221 -0.113949095 -0.0174402270 0.084352113
## Gapdh 0.07571862 0.00168726 -0.066102735 -0.0014192182 -0.005445161
## Slfn5 0.07523311 -0.03439980 0.121948327 -0.0986415897 0.015533179
## PC_6 PC_7 PC_8
## Cst7 -0.044944616 0.006920554 0.037016038
## Cd72 0.025995912 -0.028995407 0.075412988
## B2m -0.117409096 0.028275984 -0.024816009
## Ch25h 0.091269216 -0.029847492 0.101975365
## Lgals3bp -0.090008619 0.054464751 -0.040151864
## Cd52 -0.070996534 -0.009698901 0.013347348
## Fth1 -0.088306196 -0.045837095 0.002304896
## Ctsb -0.057557479 0.145524810 -0.020520516
## Ctsz -0.059787151 0.004880728 -0.053875683
## Ccl3 0.092608183 0.013953343 0.110391444
## Cd9 -0.034012305 -0.023169678 -0.038052392
## Bst2 -0.044777599 -0.015277668 -0.036204519
## Fxyd5 0.020256812 -0.047032567 0.025598323
## Ctsd -0.073937168 0.024783903 -0.049834277
## H2-K1 -0.098054435 0.023869165 -0.004797431
## Mif -0.006587504 -0.031618471 -0.033623015
## Lpl 0.060525346 0.027615059 0.042608584
## H2-D1 -0.119351346 0.042041655 -0.001854407
## Fabp5 0.028043129 0.026398805 0.049128937
## Spp1-EGFP 0.070929233 -0.050363638 0.079791275
## Lilrb4a 0.088976630 0.015712180 0.104512794
## Csf1 0.076990407 0.029524221 0.057088988
## Tspo -0.038027118 -0.051554132 -0.019477712
## Ccl2 0.085629727 -0.054186904 0.111073502
## Ccl4 0.127044371 -0.011092153 0.146989600
## Cd63 -0.057374763 0.093133316 -0.008161522
## Lyz2 -0.120994907 0.120351793 -0.018108488
## Zbp1 0.020027446 -0.046656939 0.018233516
## Oasl2 0.001815700 -0.003693521 -0.022303578
## Fn1 0.073433024 -0.049732703 0.067528229
## Irf7 -0.001806470 -0.035682258 0.009648752
## Apoe -0.110169908 0.213888729 0.111025894
## H2-Q7 -0.069706644 0.005301151 0.009538434
## Rsad2 0.078424509 -0.042403429 0.053046329
## Spp1 0.054218209 -0.018373564 0.085418704
## Cxcl16 0.006855720 -0.012366960 0.035189155
## Phf11b 0.025043245 -0.072215760 0.010316775
## Ftl1 -0.036603625 0.027822457 -0.030990694
## Gapdh -0.011628249 -0.042104589 -0.070552164
## Slfn5 0.025909198 0.014281558 -0.039929855(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(PC_1 ) %>% head(40)## PC_1 PC_2 PC_3 PC_4 PC_5
## Crybb1 -0.06463817 0.0294974727 -0.0354208164 -0.097408764 0.006862056
## Ccr5 -0.06204224 0.0110291106 0.0140474164 -0.043660322 -0.028842188
## Sox4 -0.04764682 0.0121422496 -0.0171804830 -0.006418549 -0.041590884
## Fchsd2 -0.04629952 0.0227555136 -0.0045365074 -0.049839575 -0.046466318
## Hpgd -0.03950759 -0.0059255018 0.0374309333 0.088391440 -0.037506101
## Clec4a3 -0.03504186 -0.0010099934 0.0138714112 -0.002062804 0.005103966
## Bank1 -0.03184373 0.0124236882 -0.0018493449 -0.016879860 -0.037725693
## Fcgrt -0.03149697 0.0198491582 -0.0373920652 -0.090065353 0.011228072
## Lst1 -0.03141994 0.0077358960 -0.0143684268 -0.070862561 0.049164853
## Klf7 -0.03067361 0.0135257120 -0.0053927050 -0.030175833 -0.019121449
## Mrc1 -0.03063504 0.0457878200 -0.0460962423 -0.119637080 0.036337511
## Igfbp4 -0.02825798 0.0380816918 -0.0362660542 -0.104365500 0.021842448
## Serpinf1 -0.02810006 0.0142248546 -0.0151006007 -0.037692334 0.009542200
## Camk2n1 -0.02488802 0.0103614022 0.0006635641 -0.010470197 -0.026282061
## Rgs7bp -0.02397541 0.0188955579 -0.0266946601 -0.039995923 -0.006532197
## Fry -0.02292624 0.0080681519 0.0069077812 -0.018805924 -0.024590218
## Gpr165 -0.02218366 0.0192150939 0.0041709348 -0.013638851 -0.015836455
## Cbfa2t3 -0.02218013 0.0135865919 0.0076007074 0.010575410 -0.017679377
## Tnfrsf17 -0.02130499 0.0002597742 0.0072914149 0.016976757 -0.006562451
## Il7r -0.02108119 -0.0049129162 0.0423682248 0.066751928 -0.044075832
## Cryba4 -0.01920738 0.0131492066 -0.0383922518 -0.085652922 0.041492931
## Gbp7 -0.01903686 0.0122138706 0.0314363471 -0.086906707 -0.003969222
## Filip1l -0.01816183 0.0084317125 0.0288659311 0.009828060 -0.037359629
## Lrba -0.01709218 -0.0079392318 0.0445312204 0.087127738 -0.062444147
## Ppm1l -0.01685085 -0.0006411684 0.0120430166 0.021498073 -0.027764018
## Rdx -0.01563234 0.0322523197 -0.0101063926 -0.028082107 -0.029036583
## Zbtb20 -0.01536199 0.0002726355 0.0438005938 0.048369792 -0.036211015
## Ccr1 -0.01517874 0.0077022669 -0.0120585311 -0.041937733 0.019017410
## Ddit4 -0.01486896 -0.0009223939 0.0121797329 0.023814395 -0.020442643
## Nav2 -0.01409866 -0.0142929110 0.0593419775 0.096965281 -0.061809571
## Arid4a -0.01399211 0.0061103712 0.0168647372 -0.002183105 -0.013359338
## Igf1r -0.01345657 0.0071848689 -0.0088389407 -0.006032026 -0.032765500
## Fam102b -0.01311604 -0.0077763073 0.0377571497 0.088769528 -0.083473741
## Tet1 -0.01295294 0.0037873297 0.0004408433 -0.002861765 -0.013703246
## Khdrbs3 -0.01269506 0.0065244511 0.0181147008 0.048626673 -0.051940118
## Tmem52 -0.01262413 -0.0011493721 0.0012116385 -0.023795250 0.014379272
## Gas6 -0.01261235 0.0120974723 -0.0605701313 -0.088664181 -0.006754880
## Pecam1 -0.01244409 0.0009211247 0.0080679071 0.008475348 -0.029843046
## Upk1b -0.01235602 -0.0033641098 0.0245727728 0.052303857 -0.031521103
## Tmem176a -0.01177680 0.0057389080 0.0020905264 -0.013429253 0.033399460
## PC_6 PC_7 PC_8
## Crybb1 0.0238769004 0.011048220 -0.0028997115
## Ccr5 0.0569038358 0.055592046 -0.0139333431
## Sox4 0.0435266299 0.002963676 -0.0356961419
## Fchsd2 0.0294158288 0.089019658 -0.0232120558
## Hpgd 0.0121485849 -0.105016017 -0.0845414809
## Clec4a3 0.0158861141 0.012877214 -0.0242974813
## Bank1 0.0407457605 -0.001907312 -0.0328508027
## Fcgrt 0.0137658557 0.041666066 0.0123113148
## Lst1 -0.0124565559 0.004216013 -0.0123912356
## Klf7 0.0520305912 0.022738785 -0.0090042175
## Mrc1 0.0533247606 0.121430451 0.0255994223
## Igfbp4 0.0439007800 0.021553984 0.0089468423
## Serpinf1 0.0317907236 -0.009130330 -0.0196456924
## Camk2n1 0.0175620237 -0.027500212 -0.0389277243
## Rgs7bp 0.0392545109 0.017083206 -0.0236065790
## Fry 0.0261683839 0.028172838 -0.0207550320
## Gpr165 0.0475979302 0.008874432 0.0004577951
## Cbfa2t3 0.0264010781 -0.009471263 -0.0377540937
## Tnfrsf17 0.0142571530 -0.023715601 -0.0065563969
## Il7r 0.0023577848 -0.067293856 -0.0658356935
## Cryba4 -0.0142049472 0.055312942 -0.0181824270
## Gbp7 0.0514019814 0.019025147 -0.0137293020
## Filip1l 0.0351766859 -0.010086337 -0.0313252475
## Lrba 0.0202761343 -0.079276963 -0.0676261588
## Ppm1l 0.0192811498 -0.008495790 -0.0228033747
## Rdx 0.0118061661 0.013140204 -0.0256400848
## Zbtb20 0.0048243371 0.022296347 -0.0304752079
## Ccr1 0.0099123317 0.024381823 0.0193727422
## Ddit4 0.0137738355 -0.038357488 -0.0214061887
## Nav2 -0.0108633755 -0.027479937 -0.0837771731
## Arid4a 0.0190071055 0.036486963 -0.0422321973
## Igf1r 0.0076440834 0.047229747 -0.0039480577
## Fam102b 0.0115346031 -0.071793430 -0.1078576193
## Tet1 0.0114894028 -0.003000095 -0.0217079305
## Khdrbs3 0.0148321090 -0.073473649 -0.0534153922
## Tmem52 -0.0002480091 0.011375446 0.0103857805
## Gas6 0.0296092199 0.077575629 -0.0062290112
## Pecam1 0.0091984494 -0.009386596 -0.0231653032
## Upk1b 0.0075888593 -0.053299973 -0.0420004699
## Tmem176a -0.0207360998 -0.003450578 -0.0048027751decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.2)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 23692
## Number of edges: 592986
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7501
## Number of communities: 16
## Elapsed time: 4 secondsRun UMAP/tSNE
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 287)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 15:52:47 UMAP embedding parameters a = 0.9922 b = 1.112## 15:52:47 Read 23692 rows and found 20 numeric columns## 15:52:47 Using Annoy for neighbor search, n_neighbors = 20## 15:52:47 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:52:50 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpIHLR9L\filedccc2e1a3c91
## 15:52:50 Searching Annoy index using 1 thread, search_k = 2000
## 15:52:56 Annoy recall = 100%
## 15:52:57 Commencing smooth kNN distance calibration using 1 thread
## 15:52:58 Initializing from normalized Laplacian + noise
## 15:52:59 Commencing optimization for 200 epochs, with 683800 positive edges
## 15:53:22 Optimization finished#saveRDS(GEX.seur,"GEX0811.seur.pure_test.rds")DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "preAnno") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno2")
GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(0,1,9,10,6,8,11,14,
2,3,4,7,12,13,5,15))VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0.05, group.by = "FB.new", cols = color.FBnew)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"),
ncol = 2,
pt.size = 0, group.by = "FB.new", cols = color.FBnew)
markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
"Cd74","Lyz2","Ccl4",
"Aif1","P2ry12","C1qa","Spp1",
"Top2a","Pcna","Mki67","Mcm6",
"Cx3cr1","Il4ra","Il13ra1","Spp1-EGFP",
"Fabp5","Hmox1","Ms4a7","Cenpa")FeaturePlot(GEX.seur,
features = markers.mig,
ncol = 4)
DotPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
"Fcer1g","Il4ra","Il13ra1"),
group.by = "sort_clusters")
DotPlot(GEX.seur, features = markers.mig, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
preAnno
# Anno1
GEX.seur$Anno1 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(0)] <- "MIG.a1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(1)] <- "MIG.a2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(9)] <- "MIG.a3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(10)] <- "MIG.a4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(6)] <- "MIG.a5"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(8)] <- "MIG.a6"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(11)] <- "MIG.a7"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(14)] <- "MIG.a8"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(2)] <- "MIG.b1"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(3)] <- "MIG.b2"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(4)] <- "MIG.b3"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(7)] <- "MIG.b4"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(12)] <- "MIG.b5"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(13)] <- "MIG.b6"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(5)] <- "MIG.b7"
GEX.seur$Anno1[GEX.seur$Anno1 %in% c(15)] <- "MIG.b8"
GEX.seur$Anno1 <- factor(GEX.seur$Anno1,
levels = c(paste0("MIG.a",1:8),
paste0("MIG.b",1:8)))
# Anno2
GEX.seur$Anno2 <- as.character(GEX.seur$Anno1)
GEX.seur$Anno2 <- gsub("1|2|3|4|5|6|7|8","",GEX.seur$Anno2)
GEX.seur$Anno2 <- factor(GEX.seur$Anno2,
levels = c("MIG.a","MIG.b"))pp.new.1a <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pp.new.1a
pp.new.1b <- DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 3.45, group.by = "Anno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "Anno2", cols = color.cnt[c(3,6)])
pp.new.1b
pp.new.1c <- DimPlot(GEX.seur, reduction = "umap", group.by = "FB.new", split.by = "FB.new", cols = color.FBnew, ncol = 5)
pp.new.1c 
pp.new.1d <- DimPlot(GEX.seur, reduction = "umap", group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)
pp.new.1d 
markers
# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "Anno1"
#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE,
# min.pct = 0.05,
# test.use = "MAST",
# logfc.threshold = 0.25)
GEX.markers.pre <- read.table("sc10x_LYN.0811.marker.Anno1.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)## # A tibble: 128 x 7
## # Groups: cluster [16]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 0 0.546 0.994 0.962 0 MIG.a1 Rps19
## 2 0 0.526 0.997 0.988 0 MIG.a1 Rps11
## 3 0 0.524 0.999 0.993 0 MIG.a1 Rps27a
## 4 3.01e-318 0.519 0.994 0.97 5.53e-314 MIG.a1 Rps18
## 5 1.52e-264 0.558 0.962 0.89 2.79e-260 MIG.a1 Ltc4s
## 6 5.36e-168 0.591 0.812 0.675 9.85e-164 MIG.a1 Crybb1
## 7 4.80e-155 0.566 0.657 0.516 8.81e-151 MIG.a1 2410006H16Rik
## 8 2.39e- 95 0.587 0.44 0.29 4.39e- 91 MIG.a1 Mrc1
## 9 1.18e-273 0.718 0.891 0.706 2.17e-269 MIG.a2 Ccr5
## 10 6.88e-264 0.720 0.795 0.527 1.26e-259 MIG.a2 Fchsd2
## # ... with 118 more rowsGEX.markers.anno <- GEX.markers.pre
GEX.markers.anno$cluster <- factor(as.character(GEX.markers.anno$cluster),
levels = levels(GEX.seur$Anno1))
markers.pre_t48 <- (GEX.markers.anno %>% group_by(cluster) %>%
filter(pct.1>0.1) %>%
top_n(n = 48, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t60 <- (GEX.markers.anno %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl|Gm|Rik$",GEX.markers.anno$gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t120 <- (GEX.markers.anno %>% group_by(cluster) %>%
filter(pct.1>0.1 & gene %in% grep("Rps|Rpl|Gm|Rik$",GEX.markers.anno$gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$geneDotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[1:64])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[65:128])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[129:192])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[193:256])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[257:320])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[321:384])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[385:448])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[449:512])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
DotPlot(GEX.seur, group.by = "Anno1", features = rev(markers.pre_t60[513:568])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
cell composition
usual
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.new,
clusters=GEX.seur$Anno1),
main = "Cell Count",
gaps_row = c(1,3,5,6,8),
gaps_col = c(8),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.new,
clusters=GEX.seur$Anno1)),
main = "Cell Ratio",
gaps_row = c(1,3,5,6,8),
gaps_col = c(8),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent = T)$gtable,
ncol = 1)
glm.nb-anova.LRT
levels(GEX.seur$FB.new)## [1] "P05.CTR.1" "P05.MIG.1" "P05.MIG.2" "P05.GFP.1" "P05.GFP.2" "P28.CTR.1"
## [7] "P28.MIG.1" "P28.MIG.2" "P28.GFP.1" "P28.GFP.2"GEX.seur$rep <- gsub("P05.CTR.|P05.MIG.|P05.GFP.|P28.CTR.|P28.MIG.|P28.GFP.","rep",as.character(GEX.seur$FB.new))
head(GEX.seur$rep)## AAACCCAAGAAATCCA-1 AAACCCAAGCGATCGA-1 AAACCCAAGGTAGCCA-1 AAACCCAAGGTTGAGC-1
## "rep1" "rep2" "rep1" "rep1"
## AAACCCAAGTCACTCA-1 AAACCCAAGTGGAATT-1
## "rep1" "rep2"stat_sort <- GEX.seur@meta.data[,c("cnt","rep","Anno1")]dim(stat_sort )## [1] 23692 3head(stat_sort )## cnt rep Anno1
## AAACCCAAGAAATCCA-1 P05.GFP rep1 MIG.a2
## AAACCCAAGCGATCGA-1 P05.MIG rep2 MIG.a2
## AAACCCAAGGTAGCCA-1 P28.GFP rep1 MIG.b7
## AAACCCAAGGTTGAGC-1 P28.GFP rep1 MIG.b3
## AAACCCAAGTCACTCA-1 P28.GFP rep1 MIG.b4
## AAACCCAAGTGGAATT-1 P05.MIG rep2 MIG.a4stat_sort.s <- stat_sort %>%
group_by(cnt,rep,Anno1) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup()## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.dim(stat_sort.s)## [1] 151 5#stat_sort.s#write.csv(stat_sort.s, "figures0921/2.pure/composition/composition_Anno1.stat.csv")#write.csv(stat_sort.s %>% reshape2::recast(cnt + rep ~ Anno1, measure.var = 4),
# "figures0921/2.pure/composition/composition_Anno1.stat_df.count.csv")
#write.csv(stat_sort.s %>% reshape2::recast(cnt + rep ~ Anno1, measure.var = 5),
# "figures0921/2.pure/composition/composition_Anno1.stat_df.ratio.csv")plot
pcom.sort <- stat_sort.s %>% #filter(stage %in% c("P00","P07","P14","P21","P28")) %>%
ggplot(aes(x = Anno1, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.cnt, name="") +
labs(title="Cell Composition", y="Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=Anno1, y=100*prop, color=cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=1.75,
stroke=0.15, show.legend = F)
pcom.sort
significance
some condition no replicates, not run
test some markers
marker.fig0921 <- c("Ptprc","P2ry12","Spp1","Spp1-EGFP",
"Cx3cr1","Csf1r","Aif1","C1qa",
"Ccl4","Lpl","Fabp5","Csf1",
"Ctsf","Ccr5","Arsb","Itgb1",
"Cd81","Lag3","Cd34","H2-K1",
"Top2a","Pcna","Mki67","Mcm6"
#"Cd74","H2-Aa","Ccr2","Cd209a",
#"Mrc1","Pf4","Dab2","Lyve1",
#"Tubb3","Actl6b","Syt1","Sst"
)
pp.new.2a <- FeaturePlot(GEX.seur, features = marker.fig0921, ncol = 4,raster = T, pt.size = 2.75)
pp.new.2a
pp.new.2b <- DotPlot(GEX.seur, features = marker.fig0921, group.by = "Anno1",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
pp.new.2b
pp.new.2c <- VlnPlot(GEX.seur, features = marker.fig0921, ncol = 4, pt.size = 0, raster = T) #& geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.new.2c
pp.new.2c1 <- VlnPlot(GEX.seur, features = marker.fig0921, ncol = 4, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.new.2c1
signature score
#### 10x data, calculate signature score
#
## The code below is from Adam Hamber
## 2D scoring by Itay
get_controls <- function(counts, gene.list, verbose=F, control.genes.per.gene=10)
{
# Itay: "Such scores are inevitably correlated with cell complexity so to avoid
# that I subtract a "control" score which is generated by averaging over a control
# gene set. Control gene sets are chosen to contain 100 times more genes than the
# real gene set (analogous to averaging over 100 control sets of similar size) and
# to have the same distribution of population/bulk - based expression levels as the
# real gene set, such that they are expected to have the same number of "zeros" and
# to eliminate the correlation with complexity."
# ---------------------------------------------------------------------------------
# Going to find control points by finding the closest genes in terms of expression level and % of the time we observe it
if(verbose){cat(sprintf("Finding %s background genes based on similarity to given gene set [%s genes] \n",
control.genes.per.gene*length(gene.list), length(gene.list)))}
cat("Summarizing data \n")
summary = data.frame(gene=row.names(counts), mean.expr = Matrix::rowMeans(counts), fract.zero = Matrix::rowMeans(counts==0), stringsAsFactors = F)
#summary = data.frame(gene=row.names(counts), mean.expr = apply(counts,1,mean), fract.zero = apply(counts==0,1,mean), stringsAsFactors = F)
summary$mean.expr.s = scale(summary$mean.expr)
summary$fract.zero.s = scale(summary$fract.zero)
actual.genes = summary[summary$gene %in% gene.list,]
background.genes = summary[!summary$gene %in% gene.list,]
#find the 10 closest genes to each cell cycle marker gene and add them to the lists of control genes
get_closest_genes <- function(i)
{
background.genes$dist = sqrt((background.genes$mean.expr.s - actual.genes$mean.expr.s[i])^2 +
(background.genes$fract.zero.s - actual.genes$fract.zero.s[i])^2)
ordered = background.genes$gene[order(background.genes$dist)]
ordered = ordered[!ordered %in% controls] # don't take genes that already appear in the list
closest = head(ordered, n=control.genes.per.gene)
return(closest)
}
controls = c();
for (i in 1:length(gene.list)){
#info(sprintf("Finding %s control genes for %s", control.genes.per.gene, gene.list[i]))
closest = get_closest_genes(i)
#info(sprintf("Found %s: ", length(closest)))
controls = unique(c(controls, closest))
}
if(verbose){cat(sprintf("Control gene selection complete. %s genes found. \n", length(controls)))}
#print(controls)
return(controls)
}
## Define calculate function
calculate_signature_score <- function(count_matrix, gene_list){
control_gene <- get_controls(counts = count_matrix,
gene.list = gene_list)
signature_score <- colMeans(count_matrix[gene_list, ], na.rm = TRUE) -
colMeans(count_matrix[control_gene, ], na.rm = TRUE)
return(signature_score)
}
add_geneset_score <- function(obj, geneset, setname){
score <- calculate_signature_score(as.data.frame(obj@assays[['RNA']]@data),
geneset)
obj <- AddMetaData(obj,
score,
setname)
return(obj)
}## proc_DEG()
# to process edgeR result for DEGs' comparison
# mat_cut is a matrix after advanced filtering, 'like TPM > n in at least one condition'
#
proc_DEG <- function(deg, p.cut=0.05, FC.cut = 2, padj=TRUE, abs=TRUE, mat_cut=NULL, gene_cut=NULL){
rownames(deg) <- deg$gene
if(padj==TRUE){
deg <- deg %>% filter(padj < p.cut)
}else{
deg <- deg %>% filter(pvalue < p.cut)
}
if(abs==TRUE){
deg <- deg %>% filter(abs(FC) > FC.cut)
}else if(FC.cut >0){
deg <- deg %>% filter(FC > FC.cut)
}else{
deg <- deg %>% filter(FC < FC.cut)
}
if(!is.null(mat_cut)){
deg <- deg[rownames(deg) %in% rownames(mat_cut),]
}
if(!is.null(gene_cut)){
deg <- deg[rownames(deg) %in% gene_cut,]
}
return(deg)
}SS2 DEGs
DEGs.SS2_P06 <- read.csv("../bulkDEGs/DEG_table.P06_SIMPLE.Spp1_PvsN.csv",row.names = 1)
DEGs.SS2_P06$Gene <- rownames(DEGs.SS2_P06)
#DEGs.SS2_P06DEGs.SS2_P28 <- read.csv("../bulkDEGs/DEG_table.P28_SIMPLE.Spp1_PvsN.csv",row.names = 1)
DEGs.SS2_P28$Gene <- rownames(DEGs.SS2_P28)
#DEGs.SS2_P28DEGs.SS2_P06["Spp1",]## log2FoldChange logCPM LR pvalue padj FC Gene
## Spp1 6.05263 7.836443 232.6566 1.57e-52 7.95e-49 66.37783 Spp1DEGs.SS2_P28["Spp1",]## log2FoldChange logCPM LR pvalue padj FC Gene
## Spp1 9.978018 10.00778 146.9156 8.19e-34 6.34e-32 1008.516 Spp1DEGs.list <- list(P06.Spp1_up = proc_DEG(DEGs.SS2_P06, p.cut = 0.05, FC.cut = 1.5, padj = T, abs = FALSE)$Gene,
P28.Spp1_up = proc_DEG(DEGs.SS2_P28, p.cut = 0.05, FC.cut = 1.5, padj = T, abs = FALSE)$Gene)lapply(DEGs.list, length)## $P06.Spp1_up
## [1] 177
##
## $P28.Spp1_up
## [1] 1259pp.DEG <- venn::venn(DEGs.list,
zcolor = 'style',
ggplot = T) +
labs( title=" cutoff: padj<0.05, FC>1.5")+
theme(plot.title = element_text(size=15))
pp.DEG 
DEG.comp <- list(Spp1_up.P06only = setdiff(DEGs.list$P06.Spp1_up, DEGs.list$P28.Spp1_up),
Spp1_up.both = intersect(DEGs.list$P06.Spp1_up, DEGs.list$P28.Spp1_up),
Spp1_up.P28only = setdiff(DEGs.list$P28.Spp1_up, DEGs.list$P06.Spp1_up))lapply(DEG.comp, length)## $Spp1_up.P06only
## [1] 60
##
## $Spp1_up.both
## [1] 117
##
## $Spp1_up.P28only
## [1] 1142names(DEG.comp)## [1] "Spp1_up.P06only" "Spp1_up.both" "Spp1_up.P28only"for(NN in names(DEG.comp)){
write.table(DEG.comp[[NN]],paste0("figures0921/2.pure/signature/overlap/DEGs.",NN,".num",length(DEG.comp[[NN]]),".txt"), col.names = F, row.names = F ,quote = F)
}GEX.seur <- add_geneset_score(GEX.seur, DEG.comp$Spp1_up.P06only, "Spp1_up.P06only")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DEG.comp$Spp1_up.both, "Spp1_up.both")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DEG.comp$Spp1_up.P28only, "Spp1_up.P28only")## Summarizing dataplot
samples
VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "FB.new", cols = color.FBnew, ncol = 3)
VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "cnt", cols = color.cnt, ncol = 3, pt.size = 0) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P05.MIG","P05.GFP"),
c("P28.MIG","P28.GFP")),
label.y = c(0.3,0.3),
size=3.5
)
pcomp.1a <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[1], group.by = "cnt", cols = color.cnt, ncol = 1, pt.size = 0) + NoLegend() +
labs(x="",y="Sigature Score of DEGs") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P05.MIG","P05.GFP"),
c("P05.CTR","P05.GFP"),
c("P28.MIG","P28.GFP"),
c("P28.CTR","P28.GFP")),
label.y = c(0.25,0.3,0.25,0.3),
size=2
)
pcomp.1a
pcomp.1b <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[2], group.by = "cnt", cols = color.cnt, ncol = 1, pt.size = 0) + NoLegend() +
labs(x="",y="Sigature Score of DEGs") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P05.MIG","P05.GFP"),
c("P05.CTR","P05.GFP"),
c("P28.MIG","P28.GFP"),
c("P28.CTR","P28.GFP")),
label.y = c(0.35,0.5,0.45,0.55),
size=2
)
pcomp.1b
pcomp.1c <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[3], group.by = "cnt", cols = color.cnt, ncol = 1, pt.size = 0) + NoLegend() +
labs(x="",y="Sigature Score of DEGs") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P05.MIG","P05.GFP"),
c("P05.CTR","P05.GFP"),
c("P28.MIG","P28.GFP"),
c("P28.CTR","P28.GFP")),
label.y = c(0.5,0.55,0.55,0.6),
size=2
)
pcomp.1c
#
ggsave("figures0921/2.pure/signature/sig_score.DEGs.samples.Spp1_up.P06only.pdf",
width = 3.5,height = 4,
plot = pcomp.1a)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.samples.Spp1_up.both.pdf",
width = 3.5,height = 4,
plot = pcomp.1b)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.samples.Spp1_up.P28only.pdf",
width = 3.5,height = 4,
plot = pcomp.1c)Anno1
VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "Anno1", ncol = 3)
VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), group.by = "Anno1", ncol = 3, pt.size = 0) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P05.MIG","P05.GFP"),
c("P28.MIG","P28.GFP")),
label.y = c(0.3,0.3),
size=3.5
)## Warning: Computation failed in `stat_signif()`:
## 需要TRUE/FALSE值的地方不可以用缺少值
## Warning: Computation failed in `stat_signif()`:
## 需要TRUE/FALSE值的地方不可以用缺少值
## Warning: Computation failed in `stat_signif()`:
## 需要TRUE/FALSE值的地方不可以用缺少值
pcomp.2a <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[1], group.by = "Anno1", ncol = 1, pt.size = 0) + NoLegend() +
labs(x="",y="Sigature Score of DEGs") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("MIG.a3","MIG.a5"),
c("MIG.a4","MIG.a5"),
c("MIG.a6","MIG.a5"),
c("MIG.b5","MIG.b6"),
c("MIG.b4","MIG.b6"),
c("MIG.b6","MIG.b7"),
c("MIG.b6","MIG.b8")),
label.y = c(0.25,0.2,0.17,0.25,0.32,0.23,0.28),
size=2
)
pcomp.2a
pcomp.2b <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[2], group.by = "Anno1", ncol = 1, pt.size = 0) + NoLegend() +
labs(x="",y="Sigature Score of DEGs") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("MIG.a3","MIG.a5"),
c("MIG.a4","MIG.a5"),
c("MIG.a6","MIG.a5"),
c("MIG.b5","MIG.b6"),
c("MIG.b4","MIG.b6"),
c("MIG.b6","MIG.b7"),
c("MIG.b6","MIG.b8")),
label.y = c(0.4,0.35,0.3,0.48,0.55,0.52,0.6),
size=2
)
pcomp.2b
pcomp.2c <- VlnPlot(GEX.seur, features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only")[3], group.by = "Anno1", ncol = 1, pt.size = 0) + NoLegend() +
labs(x="",y="Sigature Score of DEGs") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("MIG.a3","MIG.a5"),
c("MIG.a4","MIG.a5"),
c("MIG.a6","MIG.a5"),
c("MIG.b5","MIG.b6"),
c("MIG.b4","MIG.b6"),
c("MIG.b6","MIG.b7"),
c("MIG.b6","MIG.b8")),
label.y = c(0.55,0.5,0.47,0.575,0.6,0.55,0.59),
size=2
)
pcomp.2c
#
ggsave("figures0921/2.pure/signature/sig_score.DEGs.Anno1.Spp1_up.P06only.pdf",
width = 7.5,height = 4,
plot = pcomp.2a)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.Anno1.Spp1_up.both.pdf",
width = 7.5,height = 4,
plot = pcomp.2b)
ggsave("figures0921/2.pure/signature/sig_score.DEGs.Anno1.Spp1_up.P28only.pdf",
width = 7.5,height = 4,
plot = pcomp.2c)#saveRDS(GEX.seur,"./GEX0811.seur.pure_Anno1.rds")forVelocity
head(GEX.seur,4)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE 2998 1458 2.701801 15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE 5474 2107 5.608330 17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE 7120 2540 3.609551 12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE 4650 1705 4.430108 19.03226
## FB.info S.Score G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214 G1 2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219 G1 1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489 G1 9
## AAACCCAAGGTTGAGC-1 P28.GFP.1 0.03069876 -0.11000513 S 5
## seurat_clusters RNA_snn_res.1.2 cnt sort_clusters
## AAACCCAAGAAATCCA-1 1 1 P05.GFP 1
## AAACCCAAGCGATCGA-1 1 1 P05.MIG 1
## AAACCCAAGGTAGCCA-1 5 5 P28.GFP 5
## AAACCCAAGGTTGAGC-1 4 4 P28.GFP 4
## pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1 0.0000000 Singlet 0.0000000
## AAACCCAAGCGATCGA-1 0.4970414 Singlet 0.4911243
## AAACCCAAGGTAGCCA-1 0.4437870 Singlet 0.4260355
## AAACCCAAGGTTGAGC-1 0.1420118 Singlet 0.1183432
## DoubletFinder0.1 preAnno1 preAnno2 preAnno FB.new Anno1
## AAACCCAAGAAATCCA-1 Singlet MIG.a2 MIG.a MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1 Doublet MIG.a2 MIG.a MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1 Doublet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1 Singlet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b3
## Anno2 rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1 0.07799678 0.01572847 0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2 -0.04681886 -0.05326058 0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1 -0.03821920 0.29210607 0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1 -0.00965059 0.10195132 0.4066607output
mainly follow the official scVelo tutorial
https://smorabit.github.io/tutorials/8_velocyto/
save metadata table
## re-process cell barcodes so that next would be easier in jupyter notebook
#GEX.seur$tissue <- sapply(GEX.seur$orig.ident, function(x){base::strsplit(x,split=".",fixed = T)[[1]][2]})
#GEX.seur$barcode <- paste0(sapply(rownames(GEX.seur@meta.data), function(x){base::strsplit(x,split="-",fixed = T)[[1]][1]}),
# ".",
# GEX.seur$tissue)
GEX.seur$barcode <- rownames(GEX.seur@meta.data)
GEX.seur$UMAP_1 <- GEX.seur@reductions$umap@cell.embeddings[,1]
GEX.seur$UMAP_2 <- GEX.seur@reductions$umap@cell.embeddings[,2]
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE 2998 1458 2.701801 15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE 5474 2107 5.608330 17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE 7120 2540 3.609551 12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE 4650 1705 4.430108 19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE 2150 894 2.790698 10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE 3738 1568 7.276619 20.43874
## FB.info S.Score G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214 G1 2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219 G1 1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489 G1 9
## AAACCCAAGGTTGAGC-1 P28.GFP.1 0.03069876 -0.11000513 S 5
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644 G1 11
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337 G1 6
## seurat_clusters RNA_snn_res.1.2 cnt sort_clusters
## AAACCCAAGAAATCCA-1 1 1 P05.GFP 1
## AAACCCAAGCGATCGA-1 1 1 P05.MIG 1
## AAACCCAAGGTAGCCA-1 5 5 P28.GFP 5
## AAACCCAAGGTTGAGC-1 4 4 P28.GFP 4
## AAACCCAAGTCACTCA-1 7 7 P28.GFP 7
## AAACCCAAGTGGAATT-1 10 10 P05.MIG 10
## pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1 0.00000000 Singlet 0.00000000
## AAACCCAAGCGATCGA-1 0.49704142 Singlet 0.49112426
## AAACCCAAGGTAGCCA-1 0.44378698 Singlet 0.42603550
## AAACCCAAGGTTGAGC-1 0.14201183 Singlet 0.11834320
## AAACCCAAGTCACTCA-1 0.00000000 Singlet 0.00000000
## AAACCCAAGTGGAATT-1 0.01183432 Singlet 0.01775148
## DoubletFinder0.1 preAnno1 preAnno2 preAnno FB.new Anno1
## AAACCCAAGAAATCCA-1 Singlet MIG.a2 MIG.a MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1 Doublet MIG.a2 MIG.a MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1 Doublet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1 Singlet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b3
## AAACCCAAGTCACTCA-1 Singlet MIG.b6 MIG.b MIG P28.GFP.1 MIG.b4
## AAACCCAAGTGGAATT-1 Singlet MIG.a3 MIG.a MIG P05.MIG.2 MIG.a4
## Anno2 rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1 0.07799678 0.01572847 0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2 -0.04681886 -0.05326058 0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1 -0.03821920 0.29210607 0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1 -0.00965059 0.10195132 0.4066607
## AAACCCAAGTCACTCA-1 MIG.b rep1 0.03137777 0.29949531 0.3521602
## AAACCCAAGTGGAATT-1 MIG.a rep2 -0.02037658 -0.01007519 0.3301183
## barcode UMAP_1 UMAP_2
## AAACCCAAGAAATCCA-1 AAACCCAAGAAATCCA-1 1.1502470 4.239640
## AAACCCAAGCGATCGA-1 AAACCCAAGCGATCGA-1 -0.9110060 3.796429
## AAACCCAAGGTAGCCA-1 AAACCCAAGGTAGCCA-1 -0.1814226 -7.373271
## AAACCCAAGGTTGAGC-1 AAACCCAAGGTTGAGC-1 0.4916820 -4.882718
## AAACCCAAGTCACTCA-1 AAACCCAAGTCACTCA-1 -0.8550854 -3.813902
## AAACCCAAGTGGAATT-1 AAACCCAAGTGGAATT-1 -1.9329503 3.032112sum(duplicated(GEX.seur$barcode))## [1] 0write counts matrix
write pca
write filtered genes
figures1002
GEX.seur <- readRDS("./GEX0811.seur.pure_Anno1.rds")
GEX.seur## An object of class Seurat
## 18371 features across 23692 samples within 1 assay
## Active assay: RNA (18371 features, 1389 variable features)
## 3 dimensional reductions calculated: pca, tsne, umaphead(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE 2998 1458 2.701801 15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE 5474 2107 5.608330 17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE 7120 2540 3.609551 12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE 4650 1705 4.430108 19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE 2150 894 2.790698 10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE 3738 1568 7.276619 20.43874
## FB.info S.Score G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214 G1 2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219 G1 1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489 G1 9
## AAACCCAAGGTTGAGC-1 P28.GFP.1 0.03069876 -0.11000513 S 5
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644 G1 11
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337 G1 6
## seurat_clusters RNA_snn_res.1.2 cnt sort_clusters
## AAACCCAAGAAATCCA-1 1 1 P05.GFP 1
## AAACCCAAGCGATCGA-1 1 1 P05.MIG 1
## AAACCCAAGGTAGCCA-1 5 5 P28.GFP 5
## AAACCCAAGGTTGAGC-1 4 4 P28.GFP 4
## AAACCCAAGTCACTCA-1 7 7 P28.GFP 7
## AAACCCAAGTGGAATT-1 10 10 P05.MIG 10
## pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1 0.00000000 Singlet 0.00000000
## AAACCCAAGCGATCGA-1 0.49704142 Singlet 0.49112426
## AAACCCAAGGTAGCCA-1 0.44378698 Singlet 0.42603550
## AAACCCAAGGTTGAGC-1 0.14201183 Singlet 0.11834320
## AAACCCAAGTCACTCA-1 0.00000000 Singlet 0.00000000
## AAACCCAAGTGGAATT-1 0.01183432 Singlet 0.01775148
## DoubletFinder0.1 preAnno1 preAnno2 preAnno FB.new Anno1
## AAACCCAAGAAATCCA-1 Singlet MIG.a2 MIG.a MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1 Doublet MIG.a2 MIG.a MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1 Doublet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1 Singlet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b3
## AAACCCAAGTCACTCA-1 Singlet MIG.b6 MIG.b MIG P28.GFP.1 MIG.b4
## AAACCCAAGTGGAATT-1 Singlet MIG.a3 MIG.a MIG P05.MIG.2 MIG.a4
## Anno2 rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1 0.07799678 0.01572847 0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2 -0.04681886 -0.05326058 0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1 -0.03821920 0.29210607 0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1 -0.00965059 0.10195132 0.4066607
## AAACCCAAGTCACTCA-1 MIG.b rep1 0.03137777 0.29949531 0.3521602
## AAACCCAAGTGGAATT-1 MIG.a rep2 -0.02037658 -0.01007519 0.3301183
## DAM.sig_top50 DAM.sig_top100 DAM.sig_top250 DAM.sig_top500
## AAACCCAAGAAATCCA-1 -0.15072706 -0.18592966 -0.13817389 0.02798874
## AAACCCAAGCGATCGA-1 -0.23592021 -0.16868345 -0.06001030 0.12662573
## AAACCCAAGGTAGCCA-1 0.79968333 0.61875031 0.48767439 0.50833408
## AAACCCAAGGTTGAGC-1 0.25948931 0.20225591 0.17024206 0.31339640
## AAACCCAAGTCACTCA-1 0.90775071 0.64612166 0.38459858 0.38874980
## AAACCCAAGTGGAATT-1 -0.02039969 -0.04251335 -0.02195155 0.17200874
## Cycling_score
## AAACCCAAGAAATCCA-1 -0.072838959
## AAACCCAAGCGATCGA-1 -0.045974881
## AAACCCAAGGTAGCCA-1 -0.006387924
## AAACCCAAGGTTGAGC-1 -0.010882599
## AAACCCAAGTCACTCA-1 -0.036946642
## AAACCCAAGTGGAATT-1 -0.007985408ppnew.1
Signature plot as feature plot
ppnew.1a <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
raster = T, pt.size = 3.5,
min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all")
ppnew.1a
ppnew.1b <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
cols = c("lightgrey","red"),raster = T, pt.size = 3.5,
min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all")
ppnew.1b
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)ppnew.1c <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
#cols = c("lightgrey","red"),
#cols = rev(mapal),
raster = T, pt.size = 3.5,
min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all") & scale_color_gradientn(colors = rev(mapal), limits= c(-0.16,0.65), breaks = c(0.0,0.2,0.4,0.6))
ppnew.1c
ppnew.1d <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only","Spp1_up.both","Spp1_up.P28only"), ncol = 3,
#cols = c("lightgrey","red"),
#cols = rev(mapal),
raster = T, pt.size = 3.5,
min.cutoff = -0.18, max.cutoff = 0.8, keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
ppnew.1d
mod
using own scale, highlight the hot-regions
ppnew.1d1 <- FeaturePlot(GEX.seur,features = c("Spp1_up.P06only"), ncol = 1,
#cols = c("lightgrey","red"),
#cols = rev(mapal),
raster = T, pt.size = 3.5,
min.cutoff = -0.08, max.cutoff = 0.38) & scale_color_gradientn(colors = rev(mapal), #limits= c(-0.08,0.38),
breaks=c(0.0,0.1,0.2,0.3))
ppnew.1d1
ppnew.1d2 <- FeaturePlot(GEX.seur,features = c("Spp1_up.both"), ncol = 1,
#cols = c("lightgrey","red"),
#cols = rev(mapal),
raster = T, pt.size = 3.5,
min.cutoff = -0.16, max.cutoff = 0.46) + scale_color_gradientn(colors = rev(mapal), #limits= c(-0.16,0.46),
breaks=c(0.0,0.2,0.4))
ppnew.1d2
ppnew.1d3a <- FeaturePlot(GEX.seur,features = c("Spp1_up.P28only"), ncol = 1,
#cols = c("lightgrey","red"),
#cols = rev(mapal),
raster = T, pt.size = 3.5,
min.cutoff = -0.1, max.cutoff = 0.54) & scale_color_gradientn(colors = rev(mapal), #limits= c(-0.1,0.64),
breaks=c(0.1,0.2,0.3,0.4,0.5))
ppnew.1d3a
ppnew.1d3b <- FeaturePlot(GEX.seur,features = c("Spp1_up.P28only"), ncol = 1,
#cols = c("lightgrey","red"),
#cols = rev(mapal),
raster = T, pt.size = 3.5,
min.cutoff = -0.1, max.cutoff = 0.54) & scale_color_gradientn(colors = rev(mapal), limits= c(-0.1,0.64),
breaks=c(0.0,0.2,0.4,0.6))
ppnew.1d3b
ppnew.2
DAM signature expression as feature plot (try top 50, top100, top250, top500)
DAM.sig <- read.csv("figures1002/new/ranking_of_DAM_indicator_genes.csv")
DAM.sig[1:8,]## 锘縍anking Gene Frequence Total.score
## 1 1 Spp1 11 35.98221
## 2 2 Apoe 11 31.44704
## 3 3 Ifitm3 9 20.82901
## 4 4 Ifi27l2a 9 20.44047
## 5 5 Lyz2 11 20.22463
## 6 6 Cst7 9 19.48372
## 7 7 Cd74 7 18.24720
## 8 8 Lgals3 7 18.24180DAM.list <- list(top50=DAM.sig$Gene[1:50],
top100=DAM.sig$Gene[1:100],
top250=DAM.sig$Gene[1:250],
top500=DAM.sig$Gene[1:500])DAM.list## $top50
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a" "Lyz2" "Cst7"
## [7] "Cd74" "Lgals3" "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1" "H2-Aa" "Fxyd5"
## [19] "Ccl4" "Cybb" "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3" "Axl" "Anxa5"
## [31] "Isg15" "Lgals1" "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5" "Il1b" "Crip1"
## [43] "Slfn2" "B2m" "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1"
##
## $top100
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a" "Lyz2" "Cst7"
## [7] "Cd74" "Lgals3" "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1" "H2-Aa" "Fxyd5"
## [19] "Ccl4" "Cybb" "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3" "Axl" "Anxa5"
## [31] "Isg15" "Lgals1" "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5" "Il1b" "Crip1"
## [43] "Slfn2" "B2m" "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1" "Fcgr4" "Il1rn" "Wfdc17" "Cxcl2"
## [55] "Cxcl16" "Prdx1" "Rtp4" "H2-D1" "Pkm" "Stat1"
## [61] "Anxa2" "Gpnmb" "Zbp1" "Ftl1" "Ldha" "Npc2"
## [67] "Ly6a" "Oasl2" "Gapdh" "Prdx5" "Gbp2" "Grn"
## [73] "Ifi207" "Ifitm2" "Tlr2" "Txn1" "Phf11b" "Ctsz"
## [79] "Pianp" "Cd36" "Itgax" "Fgl2" "Ccl6" "Iigp1"
## [85] "Ccl7" "H2-K1" "Pld3" "Tpi1" "Akr1a1" "Usp18"
## [91] "Rab7b" "Tmsb10" "Cd44" "Aldoa" "Hcar2" "Acod1"
## [97] "Cd300lb" "Ctsc" "Gpr65" "H2-Q7"
##
## $top250
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a" "Lyz2" "Cst7"
## [7] "Cd74" "Lgals3" "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1" "H2-Aa" "Fxyd5"
## [19] "Ccl4" "Cybb" "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3" "Axl" "Anxa5"
## [31] "Isg15" "Lgals1" "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5" "Il1b" "Crip1"
## [43] "Slfn2" "B2m" "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1" "Fcgr4" "Il1rn" "Wfdc17" "Cxcl2"
## [55] "Cxcl16" "Prdx1" "Rtp4" "H2-D1" "Pkm" "Stat1"
## [61] "Anxa2" "Gpnmb" "Zbp1" "Ftl1" "Ldha" "Npc2"
## [67] "Ly6a" "Oasl2" "Gapdh" "Prdx5" "Gbp2" "Grn"
## [73] "Ifi207" "Ifitm2" "Tlr2" "Txn1" "Phf11b" "Ctsz"
## [79] "Pianp" "Cd36" "Itgax" "Fgl2" "Ccl6" "Iigp1"
## [85] "Ccl7" "H2-K1" "Pld3" "Tpi1" "Akr1a1" "Usp18"
## [91] "Rab7b" "Tmsb10" "Cd44" "Aldoa" "Hcar2" "Acod1"
## [97] "Cd300lb" "Ctsc" "Gpr65" "H2-Q7" "Cdkn1a" "Psat1"
## [103] "Trim30a" "Cxcl14" "Srgn" "Mmp12" "Bcl2a1b" "Tap1"
## [109] "AB124611" "Xaf1" "Ly6e" "Psme1" "Ctsl" "Sp100"
## [115] "Cxcr4" "Psmb8" "AA467197" "Pgk1" "Emp3" "Csf1"
## [121] "Mcemp1" "Gusb" "Pim1" "Ctse" "Cox4i1" "Il12b"
## [127] "Msrb1" "Npm1" "Alcam" "Psme2" "Apoc2" "Bhlhe40"
## [133] "Stmn1" "Gm4951" "Pla2g7" "Plaur" "Tor3a" "Hspe1"
## [139] "Tpt1" "Ndufa1" "Flt1" "Tgfbi" "Cox6c" "Irgm1"
## [145] "Ifit2" "Uba52" "Igf2r" "Bola2" "Ank" "Tyrobp"
## [151] "Tpm4" "Ass1" "Ms4a4c" "Ifit1" "Ybx1" "Sod2"
## [157] "Cox8a" "Fam20c" "Oas1a" "Arg1" "Ms4a7" "Smpdl3a"
## [163] "Sh3pxd2b" "Fau" "Gnas" "Phf11d" "Ehd1" "Saa3"
## [169] "Cox5a" "Atox1" "Id3" "Lrpap1" "Olr1" "Sh3bgrl3"
## [175] "Sash1" "Hint1" "Il2rg" "Ctsd" "Postn" "H2-T23"
## [181] "Ucp2" "Sdc3" "Uqcrq" "Cox6a1" "Cd300lf" "Syngr1"
## [187] "Timp2" "Atp5e" "Id2" "S100a6" "Cd14" "Tubb6"
## [193] "Anp32b" "Fcgr2b" "Cd83" "Psmb9" "Bcl2a1a" "Aprt"
## [199] "Mfsd12" "Psap" "Cox6b1" "Lilr4b" "Plac8" "Glrx"
## [205] "Scimp" "Rilpl2" "Psmb6" "Gm11808" "Chmp4b" "Actr3b"
## [211] "Ly86" "Fundc2" "Ifi211" "Hif1a" "Serf2" "Dram1"
## [217] "Parp14" "Ptgs2" "Cxcl9" "Myo5a" "Gabarap" "Cyp4f18"
## [223] "Shisa5" "Lilrb4a" "Cpd" "Iqgap1" "Slfn5" "Tnfaip2"
## [229] "Nme1" "Cotl1" "Tagln2" "Mt1" "Mpeg1" "C3"
## [235] "Ube2l6" "Clic4" "Naaa" "Gas7" "Sdcbp" "Nampt"
## [241] "Ell2" "Samhd1" "Rtcb" "Eef1g" "Hmgb2" "Gng5"
## [247] "Nfil3" "Adora1" "Tpd52" "Ptger4"
##
## $top500
## [1] "Spp1" "Apoe" "Ifitm3" "Ifi27l2a"
## [5] "Lyz2" "Cst7" "Cd74" "Lgals3"
## [9] "Lpl" "Cd52" "Ccl5" "Ccl12"
## [13] "Lgals3bp" "Cd63" "H2-Ab1" "Fn1"
## [17] "H2-Aa" "Fxyd5" "Ccl4" "Cybb"
## [21] "Bst2" "Cstb" "H2-Eb1" "Fth1"
## [25] "Vim" "Tspo" "Ctsb" "Ccl3"
## [29] "Axl" "Anxa5" "Isg15" "Lgals1"
## [33] "Ccl2" "Ifi204" "Igf1" "Ch25h"
## [37] "Mif" "Cxcl10" "Plin2" "Fabp5"
## [41] "Il1b" "Crip1" "Slfn2" "B2m"
## [45] "Irf7" "Cd72" "Capg" "Ms4a6c"
## [49] "Ifit3" "Apoc1" "Fcgr4" "Il1rn"
## [53] "Wfdc17" "Cxcl2" "Cxcl16" "Prdx1"
## [57] "Rtp4" "H2-D1" "Pkm" "Stat1"
## [61] "Anxa2" "Gpnmb" "Zbp1" "Ftl1"
## [65] "Ldha" "Npc2" "Ly6a" "Oasl2"
## [69] "Gapdh" "Prdx5" "Gbp2" "Grn"
## [73] "Ifi207" "Ifitm2" "Tlr2" "Txn1"
## [77] "Phf11b" "Ctsz" "Pianp" "Cd36"
## [81] "Itgax" "Fgl2" "Ccl6" "Iigp1"
## [85] "Ccl7" "H2-K1" "Pld3" "Tpi1"
## [89] "Akr1a1" "Usp18" "Rab7b" "Tmsb10"
## [93] "Cd44" "Aldoa" "Hcar2" "Acod1"
## [97] "Cd300lb" "Ctsc" "Gpr65" "H2-Q7"
## [101] "Cdkn1a" "Psat1" "Trim30a" "Cxcl14"
## [105] "Srgn" "Mmp12" "Bcl2a1b" "Tap1"
## [109] "AB124611" "Xaf1" "Ly6e" "Psme1"
## [113] "Ctsl" "Sp100" "Cxcr4" "Psmb8"
## [117] "AA467197" "Pgk1" "Emp3" "Csf1"
## [121] "Mcemp1" "Gusb" "Pim1" "Ctse"
## [125] "Cox4i1" "Il12b" "Msrb1" "Npm1"
## [129] "Alcam" "Psme2" "Apoc2" "Bhlhe40"
## [133] "Stmn1" "Gm4951" "Pla2g7" "Plaur"
## [137] "Tor3a" "Hspe1" "Tpt1" "Ndufa1"
## [141] "Flt1" "Tgfbi" "Cox6c" "Irgm1"
## [145] "Ifit2" "Uba52" "Igf2r" "Bola2"
## [149] "Ank" "Tyrobp" "Tpm4" "Ass1"
## [153] "Ms4a4c" "Ifit1" "Ybx1" "Sod2"
## [157] "Cox8a" "Fam20c" "Oas1a" "Arg1"
## [161] "Ms4a7" "Smpdl3a" "Sh3pxd2b" "Fau"
## [165] "Gnas" "Phf11d" "Ehd1" "Saa3"
## [169] "Cox5a" "Atox1" "Id3" "Lrpap1"
## [173] "Olr1" "Sh3bgrl3" "Sash1" "Hint1"
## [177] "Il2rg" "Ctsd" "Postn" "H2-T23"
## [181] "Ucp2" "Sdc3" "Uqcrq" "Cox6a1"
## [185] "Cd300lf" "Syngr1" "Timp2" "Atp5e"
## [189] "Id2" "S100a6" "Cd14" "Tubb6"
## [193] "Anp32b" "Fcgr2b" "Cd83" "Psmb9"
## [197] "Bcl2a1a" "Aprt" "Mfsd12" "Psap"
## [201] "Cox6b1" "Lilr4b" "Plac8" "Glrx"
## [205] "Scimp" "Rilpl2" "Psmb6" "Gm11808"
## [209] "Chmp4b" "Actr3b" "Ly86" "Fundc2"
## [213] "Ifi211" "Hif1a" "Serf2" "Dram1"
## [217] "Parp14" "Ptgs2" "Cxcl9" "Myo5a"
## [221] "Gabarap" "Cyp4f18" "Shisa5" "Lilrb4a"
## [225] "Cpd" "Iqgap1" "Slfn5" "Tnfaip2"
## [229] "Nme1" "Cotl1" "Tagln2" "Mt1"
## [233] "Mpeg1" "C3" "Ube2l6" "Clic4"
## [237] "Naaa" "Gas7" "Sdcbp" "Nampt"
## [241] "Ell2" "Samhd1" "Rtcb" "Eef1g"
## [245] "Hmgb2" "Gng5" "Nfil3" "Adora1"
## [249] "Tpd52" "Ptger4" "Eif2ak2" "Areg"
## [253] "Rsad2" "Slc31a1" "Gm2000" "Cox7c"
## [257] "Tmem256" "Cox5b" "Cyba" "Il18bp"
## [261] "Selenow" "Myl6" "H2-DMa" "Rai14"
## [265] "Dab2" "Hmox1" "Tmsb4x" "Ndufa2"
## [269] "Cd68" "Ccdc86" "Atp6v1a" "Cox7a2"
## [273] "Calm1" "Uqcrh" "Socs2" "Gpi1"
## [277] "Ubl5" "Colec12" "Ifit3b" "Scpep1"
## [281] "S100a11" "F3" "Ms4a6d" "Hacd2"
## [285] "Chil3" "Tuba1c" "Rap2b" "Gp49a"
## [289] "Milr1" "Fcgr1" "Stat2" "Atp5j2"
## [293] "Uqcr10" "Ssr4" "Bcar3" "Gm49368"
## [297] "Tmem106a" "Tubb5" "Naca" "Hspa8"
## [301] "Atp5k" "Bag1" "Sec61b" "Gyg"
## [305] "Cox7b" "Ly6c2" "Top2a" "Aldh2"
## [309] "Dpp7" "Eif3k" "Cd48" "C4b"
## [313] "Pycard" "Atp5l" "Uqcr11" "Osm"
## [317] "Prelid1" "Rnf213" "Prdx4" "Arpc1b"
## [321] "Ndufv3" "Sp110" "Ndufa13" "Abca1"
## [325] "Gm1673" "Wdr89" "Sh3bgrl" "Smdt1"
## [329] "Gatm" "Gpr84" "Slamf8" "Ccrl2"
## [333] "Tomm7" "Gas8" "Ly6i" "Bcl2a1d"
## [337] "Esd" "Dhx58" "Ctsa" "Rxrg"
## [341] "Prdx6" "Ndufc1" "Polr2f" "Sdc4"
## [345] "Atp5j" "Litaf" "Atp6v0e" "Cspg4"
## [349] "Ranbp1" "Ifi35" "Fcer1g" "Calm3"
## [353] "Rhoc" "Pde4b" "Atp5g1" "Gpx3"
## [357] "Psmb1" "Gpx1" "Eef1a1" "Ly9"
## [361] "Igtp" "H2-Q6" "Herc6" "Cd9"
## [365] "Ran" "Hebp1" "Eno1" "Cdk18"
## [369] "Eef1b2" "Gbp7" "Glipr1" "Atp6v1f"
## [373] "H2-DMb1" "Btf3" "Slc25a3" "Brdt"
## [377] "Csf2ra" "Eif3f" "Hpse" "Gm14305"
## [381] "Gm14410" "H2-Q4" "Ndufb9" "Epsti1"
## [385] "Dnaaf3" "Pf4" "S100a4" "Atp5h"
## [389] "Apobec3" "Hsp90ab1" "Tnf" "Ctss"
## [393] "Gas6" "Gbp3" "Gng12" "Nceh1"
## [397] "Mgst1" "Cfl1" "Dtnbp1" "Slc2a6"
## [401] "Eif5a" "Atp5c1" "ptchd1" "Ptchd1"
## [405] "Psma7" "Lap3" "Rbm3" "Cycs"
## [409] "Cox6a2" "Abcg1" "Prr15" "Ahnak"
## [413] "Ndufb7" "Nrp1" "Lmnb1" "Ninj1"
## [417] "Mrpl33" "Arpc3" "Phlda1" "Acp5"
## [421] "Atp5g3" "C5ar1" "Arl5c" "Parp9"
## [425] "Ndufb8" "Txndc17" "Tbca" "Gm49339"
## [429] "Tma7" "Fblim1" "Eif3h" "Psme2b"
## [433] "Mrpl30" "Arpc2" "Ptma" "Rac2"
## [437] "Ctsh" "Dcstamp" "Clec4n" "Dbi"
## [441] "H2-T22" "Sem1" "Msr1" "Bax"
## [445] "S100a10" "Fabp3" "Snrpg" "Bcl2a1"
## [449] "Serpine2" "Snrpd2" "Cd180" "Pgam1"
## [453] "2310061I04Rik" "S100a1" "Serpinb6a" "Actg1"
## [457] "Uqcc2" "Tuba1b" "Neurl2" "Gcnt2"
## [461] "Smc2" "Atp5b" "Ifih1" "Nap1l1"
## [465] "Cd274" "Rgs1" "Parp12" "Serpine1"
## [469] "Nsa2" "Cebpb" "Csf2rb" "Cfb"
## [473] "Ndufa6" "Sdf2l1" "Anxa4" "Psma5"
## [477] "Nfkbiz" "Ctla2a" "Atp5o" "Ube2l3"
## [481] "Lipa" "Pfn1" "Ndufa7" "Ndufs6"
## [485] "Psmb10" "Mapkapk2" "Aif1" "Smc4"
## [489] "Itgb1bp1" "Nme2" "Pfdn5" "Rassf3"
## [493] "1810058I24Rik" "Pgls" "Clec4d" "Mrpl54"
## [497] "Tmem160" "Pomp" "Slc7a11" "F10"lapply(DAM.list, length)## $top50
## [1] 50
##
## $top100
## [1] 100
##
## $top250
## [1] 250
##
## $top500
## [1] 500GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top50, "DAM.sig_top50")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DAM.list$top100, "DAM.sig_top100")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DAM.list$top250, "DAM.sig_top250")## Summarizing dataGEX.seur <- add_geneset_score(GEX.seur, DAM.list$top500, "DAM.sig_top500")## Summarizing datappnew.2a <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
raster = T, pt.size = 3.5,
keep.scale = "all")
ppnew.2a
ppnew.2b <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
raster = T, pt.size = 3.5,
cols = c("lightgrey","red"),
keep.scale = "all")
ppnew.2b
ppnew.2c <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
raster = T, pt.size = 3.5,
keep.scale = "all") & scale_color_gradientn(colors = rev(mapal), limits= c(-0.4,1.65), breaks = c(0.0,0.5,1.0,1.5))
ppnew.2c
ppnew.2d <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
raster = T, pt.size = 3.5,
keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
ppnew.2d
comparison
each cluster, compare CTRvsGFP
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE 2998 1458 2.701801 15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE 5474 2107 5.608330 17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE 7120 2540 3.609551 12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE 4650 1705 4.430108 19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE 2150 894 2.790698 10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE 3738 1568 7.276619 20.43874
## FB.info S.Score G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214 G1 2
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219 G1 1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489 G1 9
## AAACCCAAGGTTGAGC-1 P28.GFP.1 0.03069876 -0.11000513 S 5
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644 G1 11
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337 G1 6
## seurat_clusters RNA_snn_res.1.2 cnt sort_clusters
## AAACCCAAGAAATCCA-1 1 1 P05.GFP 1
## AAACCCAAGCGATCGA-1 1 1 P05.MIG 1
## AAACCCAAGGTAGCCA-1 5 5 P28.GFP 5
## AAACCCAAGGTTGAGC-1 4 4 P28.GFP 4
## AAACCCAAGTCACTCA-1 7 7 P28.GFP 7
## AAACCCAAGTGGAATT-1 10 10 P05.MIG 10
## pANN_0.25_0.005_1267 DoubletFinder0.05 pANN_0.25_0.005_2534
## AAACCCAAGAAATCCA-1 0.00000000 Singlet 0.00000000
## AAACCCAAGCGATCGA-1 0.49704142 Singlet 0.49112426
## AAACCCAAGGTAGCCA-1 0.44378698 Singlet 0.42603550
## AAACCCAAGGTTGAGC-1 0.14201183 Singlet 0.11834320
## AAACCCAAGTCACTCA-1 0.00000000 Singlet 0.00000000
## AAACCCAAGTGGAATT-1 0.01183432 Singlet 0.01775148
## DoubletFinder0.1 preAnno1 preAnno2 preAnno FB.new Anno1
## AAACCCAAGAAATCCA-1 Singlet MIG.a2 MIG.a MIG P05.GFP.1 MIG.a2
## AAACCCAAGCGATCGA-1 Doublet MIG.a2 MIG.a MIG P05.MIG.2 MIG.a2
## AAACCCAAGGTAGCCA-1 Doublet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b7
## AAACCCAAGGTTGAGC-1 Singlet MIG.b7 MIG.b MIG P28.GFP.1 MIG.b3
## AAACCCAAGTCACTCA-1 Singlet MIG.b6 MIG.b MIG P28.GFP.1 MIG.b4
## AAACCCAAGTGGAATT-1 Singlet MIG.a3 MIG.a MIG P05.MIG.2 MIG.a4
## Anno2 rep Spp1_up.P06only Spp1_up.both Spp1_up.P28only
## AAACCCAAGAAATCCA-1 MIG.a rep1 0.07799678 0.01572847 0.2361100
## AAACCCAAGCGATCGA-1 MIG.a rep2 -0.04681886 -0.05326058 0.2947471
## AAACCCAAGGTAGCCA-1 MIG.b rep1 -0.03821920 0.29210607 0.4436883
## AAACCCAAGGTTGAGC-1 MIG.b rep1 -0.00965059 0.10195132 0.4066607
## AAACCCAAGTCACTCA-1 MIG.b rep1 0.03137777 0.29949531 0.3521602
## AAACCCAAGTGGAATT-1 MIG.a rep2 -0.02037658 -0.01007519 0.3301183
## DAM.sig_top50 DAM.sig_top100 DAM.sig_top250 DAM.sig_top500
## AAACCCAAGAAATCCA-1 -0.15072706 -0.18592966 -0.13817389 0.02798874
## AAACCCAAGCGATCGA-1 -0.23592021 -0.16868345 -0.06001030 0.12662573
## AAACCCAAGGTAGCCA-1 0.79968333 0.61875031 0.48767439 0.50833408
## AAACCCAAGGTTGAGC-1 0.25948931 0.20225591 0.17024206 0.31339640
## AAACCCAAGTCACTCA-1 0.90775071 0.64612166 0.38459858 0.38874980
## AAACCCAAGTGGAATT-1 -0.02039969 -0.04251335 -0.02195155 0.17200874
## Cycling_score
## AAACCCAAGAAATCCA-1 -0.072838959
## AAACCCAAGCGATCGA-1 -0.045974881
## AAACCCAAGGTAGCCA-1 -0.006387924
## AAACCCAAGGTTGAGC-1 -0.010882599
## AAACCCAAGTCACTCA-1 -0.036946642
## AAACCCAAGTGGAATT-1 -0.007985408ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1
ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
same.y.lims = T,
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1
ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1
ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
same.y.lims = T,
group.by = "cnt", cols = color.cnt,
ncol = 2, pt.size = 0, raster = F) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("P05.CTR","P05.GFP"),
c("P28.CTR","P28.GFP")),
label.y = c(1,1.4),
size=3
)
ppnew.2.v1
ppnew.2.v2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
split.by = "cnt", cols = color.cnt,
ncol = 1, pt.size = 0, raster = F) #& geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v2
subset
subset ‘MIG.a - P05’ and ‘MIG.b - P28’ to check ‘CTR vs GFP’ in each cluster
GEX.P05_MIGa <- subset(GEX.seur, subset = Anno2 == "MIG.a" & cnt %in% c("P05.CTR","P05.MIG","P05.GFP"))
GEX.P05_MIGa## An object of class Seurat
## 18371 features across 11960 samples within 1 assay
## Active assay: RNA (18371 features, 1389 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapGEX.P28_MIGb <- subset(GEX.seur, subset = Anno2 == "MIG.b" & cnt %in% c("P28.CTR","P28.MIG","P28.GFP"))
GEX.P28_MIGb## An object of class Seurat
## 18371 features across 11099 samples within 1 assay
## Active assay: RNA (18371 features, 1389 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapcnta
sort(as.character(unique(GEX.P05_MIGa$Anno1)))## [1] "MIG.a1" "MIG.a2" "MIG.a3" "MIG.a4" "MIG.a5" "MIG.a6" "MIG.a7" "MIG.a8"sort(as.character(unique(GEX.P05_MIGa$cnt)))## [1] "P05.CTR" "P05.GFP" "P05.MIG"cnta <- sapply(sort(as.character(unique(GEX.P05_MIGa$Anno1))),function(x){
paste0(x,"_",sort(as.character(unique(GEX.P05_MIGa$cnt)))[c(1,3,2)])
})
cnta <- as.vector(unlist(cnta))
cnta## [1] "MIG.a1_P05.CTR" "MIG.a1_P05.MIG" "MIG.a1_P05.GFP" "MIG.a2_P05.CTR"
## [5] "MIG.a2_P05.MIG" "MIG.a2_P05.GFP" "MIG.a3_P05.CTR" "MIG.a3_P05.MIG"
## [9] "MIG.a3_P05.GFP" "MIG.a4_P05.CTR" "MIG.a4_P05.MIG" "MIG.a4_P05.GFP"
## [13] "MIG.a5_P05.CTR" "MIG.a5_P05.MIG" "MIG.a5_P05.GFP" "MIG.a6_P05.CTR"
## [17] "MIG.a6_P05.MIG" "MIG.a6_P05.GFP" "MIG.a7_P05.CTR" "MIG.a7_P05.MIG"
## [21] "MIG.a7_P05.GFP" "MIG.a8_P05.CTR" "MIG.a8_P05.MIG" "MIG.a8_P05.GFP"ctta <- list()
for(i in 1:8){
ctta[[i]] <- cnta[c(i*3-2,i*3)]
}
ctta## [[1]]
## [1] "MIG.a1_P05.CTR" "MIG.a1_P05.GFP"
##
## [[2]]
## [1] "MIG.a2_P05.CTR" "MIG.a2_P05.GFP"
##
## [[3]]
## [1] "MIG.a3_P05.CTR" "MIG.a3_P05.GFP"
##
## [[4]]
## [1] "MIG.a4_P05.CTR" "MIG.a4_P05.GFP"
##
## [[5]]
## [1] "MIG.a5_P05.CTR" "MIG.a5_P05.GFP"
##
## [[6]]
## [1] "MIG.a6_P05.CTR" "MIG.a6_P05.GFP"
##
## [[7]]
## [1] "MIG.a7_P05.CTR" "MIG.a7_P05.GFP"
##
## [[8]]
## [1] "MIG.a8_P05.CTR" "MIG.a8_P05.GFP"GEX.P05_MIGa$Anno1_cnt <- factor(paste0(as.character(GEX.P05_MIGa$Anno1),
"_",
as.character(GEX.P05_MIGa$cnt)),
levels = cnta)cntb
sort(as.character(unique(GEX.P28_MIGb$Anno1)))## [1] "MIG.b1" "MIG.b2" "MIG.b3" "MIG.b4" "MIG.b5" "MIG.b6" "MIG.b7" "MIG.b8"sort(as.character(unique(GEX.P28_MIGb$cnt)))## [1] "P28.CTR" "P28.GFP" "P28.MIG"cntb <- sapply(sort(as.character(unique(GEX.P28_MIGb$Anno1))),function(x){
paste0(x,"_",sort(as.character(unique(GEX.P28_MIGb$cnt)))[c(1,3,2)])
})
cntb <- as.vector(unlist(cntb))
cntb## [1] "MIG.b1_P28.CTR" "MIG.b1_P28.MIG" "MIG.b1_P28.GFP" "MIG.b2_P28.CTR"
## [5] "MIG.b2_P28.MIG" "MIG.b2_P28.GFP" "MIG.b3_P28.CTR" "MIG.b3_P28.MIG"
## [9] "MIG.b3_P28.GFP" "MIG.b4_P28.CTR" "MIG.b4_P28.MIG" "MIG.b4_P28.GFP"
## [13] "MIG.b5_P28.CTR" "MIG.b5_P28.MIG" "MIG.b5_P28.GFP" "MIG.b6_P28.CTR"
## [17] "MIG.b6_P28.MIG" "MIG.b6_P28.GFP" "MIG.b7_P28.CTR" "MIG.b7_P28.MIG"
## [21] "MIG.b7_P28.GFP" "MIG.b8_P28.CTR" "MIG.b8_P28.MIG" "MIG.b8_P28.GFP"cttb <- list()
for(i in 1:8){
cttb[[i]] <- cntb[c(i*3-2,i*3)]
}
cttb## [[1]]
## [1] "MIG.b1_P28.CTR" "MIG.b1_P28.GFP"
##
## [[2]]
## [1] "MIG.b2_P28.CTR" "MIG.b2_P28.GFP"
##
## [[3]]
## [1] "MIG.b3_P28.CTR" "MIG.b3_P28.GFP"
##
## [[4]]
## [1] "MIG.b4_P28.CTR" "MIG.b4_P28.GFP"
##
## [[5]]
## [1] "MIG.b5_P28.CTR" "MIG.b5_P28.GFP"
##
## [[6]]
## [1] "MIG.b6_P28.CTR" "MIG.b6_P28.GFP"
##
## [[7]]
## [1] "MIG.b7_P28.CTR" "MIG.b7_P28.GFP"
##
## [[8]]
## [1] "MIG.b8_P28.CTR" "MIG.b8_P28.GFP"GEX.P28_MIGb$Anno1_cnt <- factor(paste0(as.character(GEX.P28_MIGb$Anno1),
"_",
as.character(GEX.P28_MIGb$cnt)),
levels = cntb)newplot
ppnew.2.v2a <- VlnPlot(GEX.P05_MIGa, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
same.y.lims = T,
group.by = "Anno1_cnt", cols = rep(color.cnt[1:3],8),
ncol = 2, pt.size = 0, raster = F) &
geom_boxplot(outlier.size = 0, fill="white", width=0.35, size=0.15, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.45) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = ctta,
label.y = rep(0.53,8),
size=2.1
) & theme(axis.text.x = element_text(size=7.5)) & ylim(c(-0.4,0.6))
ppnew.2.v2a
ppnew.2.v2b <- VlnPlot(GEX.P28_MIGb, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"),
same.y.lims = T,
group.by = "Anno1_cnt", cols = c(rep(color.cnt[c(4,6,5)],7),color.cnt[6:5]),
ncol = 2, pt.size = 0, raster = F) &
geom_boxplot(outlier.size = 0, fill="white", width=0.35, size=0.15, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.45) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = cttb[c(1:5)],
label.y = rep(1.55,8),
size=2.1
) & theme(axis.text.x = element_text(size=7.5)) & ylim(c(-0.5,1.75))
ppnew.2.v2b
ppnew.3
New dotplot (gene list attached)
new.marker <- as.vector(unlist(read.csv("figures1002/new/dotplot_gene.csv")))
new.marker <- unique(new.marker)
new.marker## [1] "Tmem119" "P2ry12" "P2ry13" "Cx3cr1" "Crybb1" "Sall1" "Mafb"
## [8] "Igf1" "Apoe" "Cd63" "Abcg1" "Lpl" "Psat1" "Csf1"
## [15] "Lilrb4a" "Ftl1" "Mt1" "Abca1" "Gas6" "Cxcl10" "Ccl12"
## [22] "Ccl2" "Oasl2" "Irf7" "Oasl1" "Tnf" "Il1b" "Bcl2a1b"
## [29] "Bcl2a1d" "Bcl2a1a" "Gpx4" "Bcl2" "Ccl5"## mod color
# scale_color_gradientn(colours = material.heat(100))
# Cell2020
material.heat = function(n){
colorRampPalette(
c(
"#283593", # indigo 800
"#3F51B5", # indigo
"#2196F3", # blue
"#00BCD4", # cyan
"#4CAF50", # green
"#8BC34A", # light green
"#CDDC39", # lime
"#FFEB3B", # yellow
"#FFC107", # amber
"#FF9800", # orange
"#FF5722" # deep orange
)
)(n)
}
# Immunity2019, na gray
colors.Immunity <-c("#191970","#121285","#0C0C9A","#0707B0","#0101C5","#0014CF","#0033D3","#0053D8","#0072DD","#0092E1","#00B2E6",
"#00D1EB","#23E8CD","#7AF17B","#D2FA29","#FFEB00","#FFC300","#FF9B00","#FF8400","#FF7800","#FF6B00","#FF5F00","#FF5300",
"#FF4700","#F73B00","#EF2E00","#E62300","#DD1700","#D50B00","#CD0000")ppnew.3a1 <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1"#,cols = c("midnightblue","darkorange1")
) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3a1
ppnew.3a2 <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3a2
ppnew.3b <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
cols = c("midnightblue","darkorange1")) +
scale_color_gradientn(colours = rev(mapal)) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3b
ppnew.3c <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
cols = c("midnightblue","darkorange1")) +
scale_color_gradientn(colours = material.heat(100)) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3c
ppnew.3d <- DotPlot(GEX.seur, features = new.marker, group.by = "Anno1",
cols = c("midnightblue","darkorange1")) +
scale_color_gradientn(colours = colors.Immunity) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Anno1")
ppnew.3d
ppnew.4
Cycling score as feature plot
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
CC_gene## [1] "Mcm5" "Pcna" "Tyms" "Fen1" "Mcm7" "Mcm4"
## [7] "Rrm1" "Ung" "Gins2" "Mcm6" "Cdca7" "Dtl"
## [13] "Prim1" "Uhrf1" "Cenpu" "Hells" "Rfc2" "Polr1b"
## [19] "Nasp" "Rad51ap1" "Gmnn" "Wdr76" "Slbp" "Ccne2"
## [25] "Ubr7" "Pold3" "Msh2" "Atad2" "Rad51" "Rrm2"
## [31] "Cdc45" "Cdc6" "Exo1" "Tipin" "Dscc1" "Blm"
## [37] "Casp8ap2" "Usp1" "Clspn" "Pola1" "Chaf1b" "Mrpl36"
## [43] "E2f8" "Hmgb2" "Cdk1" "Nusap1" "Ube2c" "Birc5"
## [49] "Tpx2" "Top2a" "Ndc80" "Cks2" "Nuf2" "Cks1b"
## [55] "Mki67" "Tmpo" "Cenpf" "Tacc3" "Pimreg" "Smc4"
## [61] "Ccnb2" "Ckap2l" "Ckap2" "Aurkb" "Bub1" "Kif11"
## [67] "Anp32e" "Tubb4b" "Gtse1" "Kif20b" "Hjurp" "Cdca3"
## [73] "Jpt1" "Cdc20" "Ttk" "Cdc25c" "Kif2c" "Rangap1"
## [79] "Ncapd2" "Dlgap5" "Cdca2" "Cdca8" "Ect2" "Kif23"
## [85] "Hmmr" "Aurka" "Psrc1" "Anln" "Lbr" "Ckap5"
## [91] "Cenpe" "Ctcf" "Nek2" "G2e3" "Gas2l3" "Cbx5"
## [97] "Cenpa"GEX.seur <- add_geneset_score(GEX.seur, CC_gene, "Cycling_score")## Summarizing datappnew.4a <- FeaturePlot(GEX.seur,features = c("Cycling_score"), ncol = 1,
raster = T, pt.size = 3.5)
ppnew.4a
ppnew.4b <- FeaturePlot(GEX.seur,features = c("Cycling_score"), ncol = 1,
raster = T, pt.size = 3.5,
cols = c("lightgrey","red"))
ppnew.4b
ppnew.4c <- FeaturePlot(GEX.seur,features = c("Cycling_score"), ncol = 1,
raster = T, pt.size = 3.5,
cols = c("lightgrey","red"))& scale_color_gradientn(colors = rev(mapal), breaks = c(0.00,0.25,0.50,0.75))
ppnew.4c
#saveRDS(GEX.seur,"./GEX0811.seur.pure_Anno1.rds")forGEO
GEX.seur <- readRDS("./GEX0811.seur.pure_Anno1.rds")
GEX.seur## An object of class Seurat
## 18371 features across 23692 samples within 1 assay
## Active assay: RNA (18371 features, 1389 variable features)
## 3 dimensional reductions calculated: pca, tsne, umapGEX.seur@meta.data[,grep("snn|pANN",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGAAATCCA-1 plus_SIMPLE 2998 1458 2.701801 15.01001
## AAACCCAAGCGATCGA-1 plus_SIMPLE 5474 2107 5.608330 17.97589
## AAACCCAAGGTAGCCA-1 plus_SIMPLE 7120 2540 3.609551 12.48596
## AAACCCAAGGTTGAGC-1 plus_SIMPLE 4650 1705 4.430108 19.03226
## AAACCCAAGTCACTCA-1 plus_SIMPLE 2150 894 2.790698 10.04651
## AAACCCAAGTGGAATT-1 plus_SIMPLE 3738 1568 7.276619 20.43874
## FB.info S.Score G2M.Score Phase seurat_clusters
## AAACCCAAGAAATCCA-1 P05.GFP.1 -0.06682298 -0.09221214 G1 1
## AAACCCAAGCGATCGA-1 P05.MIG.2 -0.09983915 -0.03784219 G1 1
## AAACCCAAGGTAGCCA-1 P28.GFP.1 -0.07149868 -0.02169489 G1 5
## AAACCCAAGGTTGAGC-1 P28.GFP.1 0.03069876 -0.11000513 S 4
## AAACCCAAGTCACTCA-1 P28.GFP.1 -0.01015567 -0.02782644 G1 7
## AAACCCAAGTGGAATT-1 P05.MIG.2 -0.02148368 -0.04812337 G1 10
## cnt sort_clusters DoubletFinder0.05 DoubletFinder0.1
## AAACCCAAGAAATCCA-1 P05.GFP 1 Singlet Singlet
## AAACCCAAGCGATCGA-1 P05.MIG 1 Singlet Doublet
## AAACCCAAGGTAGCCA-1 P28.GFP 5 Singlet Doublet
## AAACCCAAGGTTGAGC-1 P28.GFP 4 Singlet Singlet
## AAACCCAAGTCACTCA-1 P28.GFP 7 Singlet Singlet
## AAACCCAAGTGGAATT-1 P05.MIG 10 Singlet Singlet
## preAnno1 preAnno2 preAnno FB.new Anno1 Anno2 rep
## AAACCCAAGAAATCCA-1 MIG.a2 MIG.a MIG P05.GFP.1 MIG.a2 MIG.a rep1
## AAACCCAAGCGATCGA-1 MIG.a2 MIG.a MIG P05.MIG.2 MIG.a2 MIG.a rep2
## AAACCCAAGGTAGCCA-1 MIG.b7 MIG.b MIG P28.GFP.1 MIG.b7 MIG.b rep1
## AAACCCAAGGTTGAGC-1 MIG.b7 MIG.b MIG P28.GFP.1 MIG.b3 MIG.b rep1
## AAACCCAAGTCACTCA-1 MIG.b6 MIG.b MIG P28.GFP.1 MIG.b4 MIG.b rep1
## AAACCCAAGTGGAATT-1 MIG.a3 MIG.a MIG P05.MIG.2 MIG.a4 MIG.a rep2
## Spp1_up.P06only Spp1_up.both Spp1_up.P28only DAM.sig_top50
## AAACCCAAGAAATCCA-1 0.07799678 0.01572847 0.2361100 -0.15072706
## AAACCCAAGCGATCGA-1 -0.04681886 -0.05326058 0.2947471 -0.23592021
## AAACCCAAGGTAGCCA-1 -0.03821920 0.29210607 0.4436883 0.79968333
## AAACCCAAGGTTGAGC-1 -0.00965059 0.10195132 0.4066607 0.25948931
## AAACCCAAGTCACTCA-1 0.03137777 0.29949531 0.3521602 0.90775071
## AAACCCAAGTGGAATT-1 -0.02037658 -0.01007519 0.3301183 -0.02039969
## DAM.sig_top100 DAM.sig_top250 DAM.sig_top500 Cycling_score
## AAACCCAAGAAATCCA-1 -0.18592966 -0.13817389 0.02798874 -0.072838959
## AAACCCAAGCGATCGA-1 -0.16868345 -0.06001030 0.12662573 -0.045974881
## AAACCCAAGGTAGCCA-1 0.61875031 0.48767439 0.50833408 -0.006387924
## AAACCCAAGGTTGAGC-1 0.20225591 0.17024206 0.31339640 -0.010882599
## AAACCCAAGTCACTCA-1 0.64612166 0.38459858 0.38874980 -0.036946642
## AAACCCAAGTGGAATT-1 -0.04251335 -0.02195155 0.17200874 -0.007985408#saveRDS(GEX.seur,"I:/Shared_win/projects/202310_Spp1DAM/forGEO/seur_obj/Spp1GFP_SIMPLE.final.seur_obj.rds")